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

Last change on this file since 31066 was 31066, checked in by Ivan Raikov, 6 years ago

picnic: modifications to polynomial sampling of curves

File size: 40.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.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        (define (pointCoord axis p)
265          (coord (inexact->exact axis) p))
266
267        ;; Samples a parametric curve at regular intervals in the range xmin..xmax inclusive.
268        (define (sample-uniform)
269          (lambda (n)
270            (let* (
271                   (xmin  0.0)
272                   (xmax  1.0)
273                   (delta (- xmax xmin))
274                   (dx    (if (zero? delta) 0
275                              (if (< n 2)
276                                (error 'sample-uniform "number of iterations must be >= 2")
277                                (/ (- xmax xmin) (- n 1)))))
278                   )
279              (list-tabulate n (lambda (i) (+ xmin (* dx i))))
280              )))
281
282
283        ;; Samples a parametric curve according to a polynomial
284        ;; c0 + c1 x + c2 x^2 in the range xmin..xmax inclusive.
285        (define (sample-polynomial c0 c1 c2)
286          (lambda (n)
287            (let* (
288                   (f (lambda (x) (+ c0 (* x c1) (* x x c2))))
289                   (dx (/ 1.0 n))
290                   )
291              (let recur ((i 0.0) (fi 0.0) (lst '()))
292                (let ((fi1 (+ fi (f i))))
293                  (if (<= fi1 1.0)
294                      (recur (+ i dx) fi1 (cons fi1 lst))
295                      (reverse lst))))
296              )))
297
298
299        (define (compose-curves c1 c2)
300          (let ((c (pcurve:compose-curve (list + + +) c1 c2)))
301;            (print "compose-curves: c1 = ") (pp (pcurve:iterate-curve c1 10))
302;            (print "compose-curves: c2 = ") (pp (pcurve:iterate-curve c2 10))
303;            (print "compose-curves: c = ") (pp (pcurve:iterate-curve c 10))
304            c))
305             
306
307        (define (make-simple-curve fx fy fz n) 
308          (let ((c (pcurve:simple-curve n 1 (list fx fy fz) 0.0 1.0)))
309            c
310            ))
311       
312
313        (define (make-harmonic axis amp period phase n)
314          (let* ((freq (/ (* 2 PI) (/ 1.0 period)))
315                 (pf   (lambda (t) (* amp (cos (+ (* freq t) phase)))))
316                 (zf   (lambda (t) 0.0))
317                 (c (pcurve:simple-curve 
318                    (inexact->exact n) 1
319                    (cond ((= axis 0)
320                           (list pf zf zf))
321                          ((= axis 1)
322                           (list zf pf zf))
323                          ((= axis 2)
324                           (list pf zf zf))
325                          (else
326                           (error 'make-harmonic "invalid axis" axis))
327                          )
328                    0.0 1.0)))
329            c
330            ))
331       
332
333        (define (make-line-segment p dx dy dz) 
334          (let ((c (pcurve:line-segment 3 (list dx dy dz))))
335            (pcurve:translate-curve (point->list p) c)
336            ))
337       
338       
339        ;; auxiliary function to create segment points
340        (define (make-segment si np sp xyz)
341          (list-tabulate 
342           np
343           (lambda (i) (let* ((pi (+ sp i))
344                              (p (make-point 
345                                  (f64vector-ref (car xyz) pi)
346                                  (f64vector-ref (cadr xyz) pi)
347                                  (f64vector-ref (caddr xyz) pi))))
348                         (list si p)
349                         ))
350           ))
351
352
353         
354
355        ;;
356        ;; Creates a process of the given segments and number of points per
357        ;; segment from the given curve.
358        ;;
359        (define (make-segmented-process c f ns np)
360          (let ((xyz ((pcurve:sample-curve* c) (f (* ns np)))))
361            (let ((np1 (floor (/ (f64vector-length (car xyz)) ns))))
362              (let recur ((si ns) (ax '()))
363                (if (positive? si)
364                    (let ((sp (* (- si 1) np1)))
365                      (recur (- si 1) (append (make-segment si np1 sp xyz) ax)))
366                    ax)
367                ))
368            ))
369
370        ;;
371        ;; Non-segmented process
372        ;;
373        (define (make-process c f np)
374          (let ((xyz ((pcurve:sample-curve* c) (f np))))
375            (list-tabulate 
376             (f64vector-length (car xyz))
377             (lambda (i) 
378               (make-point 
379                (f64vector-ref (car xyz) i)
380                (f64vector-ref (cadr xyz) i)
381                (f64vector-ref (caddr xyz) i)))
382             ))
383          )
384       
385       
386        ;;
387        ;; Creates a named section containing points from the given segmented processes.
388        ;;
389        (define (make-segmented-section gid p label ps)
390          `(,label . 
391                   ,(fold (lambda (i+proc ax)
392                            (let ((seci (car i+proc)) 
393                                  (proc (cadr i+proc)))
394                              (cons
395                               `(,seci . 
396                                       ,(map (lambda (sp)
397                                               (let ((segi (car sp))
398                                                     (dpoint (cadr sp)))
399                                                 (let ((soma-distance (sqrt (dist2 p dpoint))))
400                                                   (make-genpoint dpoint gid soma-distance segi))
401                                                 ))
402                                             proc))
403                               ax)
404                              ))
405                          '() ps)
406                   ))
407       
408        ;;
409        ;; Non-segmented named section
410        ;;
411        (define (make-section gid p label ps)
412          `(,label . 
413                   ,(fold (lambda (i+proc ax)
414                            (let* ((seci (car i+proc)) 
415                                   (proc (cadr i+proc))
416                                   (pts (map (lambda (dpoint)
417                                               (let ((soma-distance (sqrt (dist2 p dpoint))))
418                                                 (make-genpoint dpoint gid soma-distance #f)))
419                                             proc)))
420                              (d "make-section: label = ~A gid = ~A p = ~A pts = ~A~%" 
421                                 label gid p proc)
422                              (cons `(,seci . ,pts) ax)
423                              ))
424                          '() ps)
425                   ))
426       
427       
428       
429        (define (make-gen-kd-tree data #!key (threshold 0.0))
430         
431          (define (update-bb p x-min y-min z-min x-max y-max z-max)
432            (let ((x (coord 0 p)) (y (coord 1 p)) (z (coord 2 p)))
433              (if (< x (x-min)) (x-min x))
434              (if (< y (y-min)) (y-min y))
435              (if (< z (z-min)) (z-min z))
436              (if (> x (x-max)) (x-max x))
437              (if (> y (y-max)) (y-max y))
438              (if (> z (z-max)) (z-max z))
439              ))
440
441
442          (let* (
443                 (t (list->kd-tree
444                     (filter (lambda (x) (<= threshold (cdr x))) data)
445                     make-value: (lambda (i v) (cdr v))
446                     make-point: (lambda (v) (apply make-point (car v)))
447                     offset: 0
448                     ))
449                 (n (kd-tree-size t))
450                 (x-min (make-parameter +inf.0))
451                 (y-min (make-parameter +inf.0))
452                 (z-min (make-parameter +inf.0))
453                 (x-max (make-parameter -inf.0))
454                 (y-max (make-parameter -inf.0))
455                 (z-max (make-parameter -inf.0))
456                 (bb (begin (kd-tree-for-each
457                             (lambda (p) (update-bb p x-min y-min z-min
458                                                    x-max y-max z-max))
459                             t)
460                            (list (x-min) (y-min) (z-min) (x-max) (y-max) (z-max))))
461                 )
462
463            (cons bb t)
464
465            ))
466
467
468
469
470
471        (define (genpoint-projection prefix my-comm myrank size cells fibers zone cell-start fiber-start)
472
473          (d "rank ~A: prefix = ~A zone = ~A~%" myrank prefix zone)
474
475          (fold (lambda (cell ax)
476
477                  (let* ((i (+ cell-start (car cell)))
478                         (root (modulo i size))
479                         (sections (cadr cell)))
480                   
481                    (fold 
482                     
483                     (lambda (sec ax)
484                       
485                       (let ((seci (car sec))
486                             (gxs  (cdr sec)))
487                         
488                         (let ((query-data
489                                (fold 
490                                 (lambda (gd ax)
491                                   (d "rank ~A: querying point ~A (~A)~%" 
492                                      myrank (genpoint-coords gd) 
493                                      (genpoint-parent-index gd))
494                                   (fold
495                                    (lambda (x ax) 
496                                      (let (
497                                            (source (car x))
498                                            (target i)
499                                            (distance (cadr x))
500                                            (segi (genpoint-segment gd))
501                                            )
502                                        (if (not segi)
503                                            (error 'genpoint-projection "point does not have segment index" gd))
504                                        (append (list source target distance segi seci) ax)
505                                        ))
506                                    ax
507                                    (delete-duplicates
508                                     (map (lambda (x) 
509                                            (d "rank ~A: query result = ~A (~A) (~A) ~%" 
510                                               myrank (kdnn-point x) (kdnn-distance x) (kdnn-parent-index x))
511                                            (list (+ fiber-start (kdnn-parent-index x))
512                                                  (+ (kdnn-distance x) (genpoint-parent-distance gd)  (kdnn-parent-distance x))
513                                                  ))
514                                          (kd-tree-near-neighbors* fibers zone (genpoint-coords gd)))
515                                     (lambda (u v)
516                                       (= (car u) (car v)))
517                                     )
518                                    ))
519                                 '() gxs)))
520                           (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
521
522                                  (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
523                             
524                             (append res1 ax))
525                           
526                           ))
527                       )
528                     ax sections)
529                    ))
530                '() cells ))
531
532
533       
534
535        (define (point-projection prefix my-comm myrank size pts fibers zone point-start nn-filter)
536          (fold (lambda (px ax)
537
538                  (d "~A: rank ~A: px = ~A zone=~A~%"  prefix myrank px zone)
539
540                  (let* ((i (+ point-start (car px)))
541                         (root (modulo i size))
542                         (dd (d "~A: rank ~A: querying point ~A (root ~A)~%" prefix myrank px root))
543                         (query-data
544                          (let ((pd (cadr px))) 
545                            (fold
546                             (lambda (x ax) 
547                               (let ((source (car x))
548                                     (target i)
549                                     (distance (cadr x)))
550                                 (if (and (> distance  0.) (not (= source target)))
551                                     (append (list source target distance) ax)
552                                     ax)
553                                 ))
554                             '()
555                             (delete-duplicates
556                              (map (lambda (x) 
557                                     (let ((res (list (car (cadar x)) 
558                                                      (+ (caddr x) (cadr (cadar x))))))
559                                       (d "~A: x = ~A res = ~A~%" prefix x res)
560                                       res))
561                                   (nn-filter pd (kd-tree-near-neighbors* fibers zone pd))
562                                   )
563                              (lambda (u v) (= (car u) (car v)))
564                              )
565                             ))
566                          ))
567
568
569                    (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
570                           (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
571                      (append res1 ax))
572                    ))
573
574                '() pts))
575       
576
577
578
579
580        (define-datatype layer-boundary layer-boundary?
581          (Bounds (top number?) (left number?) (bottom number?) (right number?))
582          (BoundsXZ (top number?) (left number?) (bottom number?) (right number?)
583                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
584          (BoundsYZ (top number?) (left number?) (bottom number?) (right number?)
585                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
586          )
587
588
589
590        (define (boundary-z-extent-function boundary)
591          (cases layer-boundary boundary
592                 (Bounds (top left bottom right) 
593                         (lambda (x y) 0.))
594                 (BoundsXZ (top left bottom right n k x y d d2) 
595                           (lambda (xp yp) 
596                             (let-values (((y0tab y1tab y2tab res)
597                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
598                               (f64vector-ref y0tab 0))))
599                 (BoundsYZ (top left bottom right n k x y d d2) 
600                           (lambda (xp yp) 
601                             (let-values (((y0tab y1tab y2tab res)
602                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
603                               (f64vector-ref y0tab 0))))
604                 ))
605
606
607        (define (point2d-rejection top left bottom right)
608            (lambda (p)
609              (let ((x (coord 0 p)) (y (coord 1 p)))
610                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
611            )
612
613
614        (define (compute-point3d p zu z-extent-function)
615          (let* ((x (coord 0 p))
616                 (y (coord 1 p))
617                 (z-extent (z-extent-function x y)))
618            (if (zero? zu)
619                (make-point x y zu)
620                (make-point x y (fp* zu z-extent)))
621            ))
622
623
624        (define (cluster-point-rejection p cluster-pts mean-distance)
625          (let ((D (* 4 mean-distance mean-distance))
626                (nn (kd-tree-nearest-neighbor cluster-pts p)))
627            (and (< (dist2 p nn) D) p)))
628
629
630
631        (define (XZAxis n k x-points z-points boundary)
632         
633          (if (not (>= n 3)) 
634              (error 'generate-boundary "XZAxis requires at least 3 interpolation points (n >= 3)"))
635         
636          (cases layer-boundary boundary
637                 (Bounds (top left bottom right) 
638
639                         (let-values (((d d2 constr errc diagn)
640                                       (bvsp-spline:compute n k x-points z-points)))
641                           
642                           (if (not (zero? errc)) 
643                               (error 'generate-boundary "error in constructing spline from boundary points" errc))
644                           
645                           (BoundsXZ top left bottom right n k x-points z-points d d2)))
646                 
647                 (else (error 'generate-boundary "boundary argument to XZAxis is already a pseudo-3D boundary")))
648          )
649
650
651        (define (Grid x-spacing y-spacing z-spacing boundary)
652
653          (match-let (((top left bottom right)
654                       (cases layer-boundary boundary
655                                   (Bounds (top left bottom right) 
656                                           (list top left bottom right))
657                                   (BoundsXZ (top left bottom right n k x y d d2) 
658                                             (list top left bottom right))
659                                   (BoundsYZ (top left bottom right n k x y d d2) 
660                                             (list top left bottom right))
661                                   )))
662
663          (let* (
664                 (x-extent   (- right left))
665                 (y-extent   (- top bottom))
666                 (z-extent-function
667                  (boundary-z-extent-function boundary))
668                 (compute-grid-points3d
669                  (lambda (p z-spacing z-extent-function)
670                    (let* ((x (coord 0 p))
671                           (y (coord 1 p))
672                           (z-extent (z-extent-function x y)))
673                      (let recur ((points '()) (z (/ z-spacing 2.)))
674                        (if (> z z-extent)
675                            points
676                            (recur (cons (make-point x y z) points) (+ z z-spacing)))
677                        ))
678                    ))
679                 )
680           
681            (d "Grid: x-spacing = ~A~%" x-spacing)
682            (d "Grid: y-spacing = ~A~%" y-spacing)
683            (d "Grid: z-spacing = ~A~%" z-spacing)
684           
685            (d "Grid: x-extent = ~A~%" x-extent)
686            (d "Grid: y-extent = ~A~%" y-extent)
687           
688            (let ((x-points (let recur ((points '()) (x (/ x-spacing 2.)))
689                              (if (> x x-extent)
690                                  (list->f64vector (reverse points))
691                                  (recur (cons x points) (+ x x-spacing)))))
692                  (y-points (let recur ((points '()) (y (/ y-spacing 2.)))
693                              (if (> y y-extent)
694                                  (list->f64vector (reverse points))
695                                  (recur (cons y points) (+ y y-spacing)))))
696                  )
697             
698              (let ((nx (f64vector-length x-points))
699                    (ny (f64vector-length y-points))
700                    )
701               
702                (let recur ((i 0) (j 0) (ax '()))
703                  (if (< i nx)
704                      (let ((x (f64vector-ref x-points i)))
705                        (if (< j ny)
706                            (let* ((y   (f64vector-ref y-points j))
707                                   (p   (make-point x y))
708                                   (p3ds (if (zero? z-spacing)
709                                             (list (make-point x y 0.0))
710                                             (compute-grid-points3d p z-spacing z-extent-function)))
711                                   )
712                              (recur i (+ 1 j) (if p3ds (append p3ds ax) ax)))
713                            (recur (+ 1 i) 0 ax)))
714                      (let* ((t (list->kd-tree ax))
715                             (n (kd-tree-size t)))
716                        (list t boundary)
717                        ))
718                  )))
719            ))
720          )
721
722        (define (UniformRandomPointProcess n x-seed y-seed boundary)
723
724          (match-let (((top left bottom right)
725                       (cases layer-boundary boundary
726                                   (Bounds (top left bottom right) 
727                                           (list top left bottom right))
728                                   (BoundsXZ (top left bottom right n k x y d d2) 
729                                             (list top left bottom right))
730                                   (BoundsYZ (top left bottom right n k x y d d2) 
731                                             (list top left bottom right))
732                                   )))
733
734          (let* (
735                 (x-extent   (- right left))
736                 (y-extent   (- top bottom))
737                 (z-extent-function (boundary-z-extent-function boundary))
738                 )
739
740            (let ((x-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init x-seed)))
741                  (y-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init y-seed)))
742                  (z-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init (current-milliseconds)))))
743             
744              (let ((point-rejection1 (point2d-rejection top left bottom right)))
745               
746                (let recur ((i 0) (ax '()))
747                  (if (< i n)
748                      (let ((x (f64vector-ref x-points i))
749                            (y (f64vector-ref y-points i))
750                            (z (f64vector-ref z-points i)))
751                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
752                          (let ((p3d 
753                                 (and (point-rejection1 p)
754                                      (compute-point3d p z z-extent-function))))
755                            (recur (+ 1 i) (if p3d (cons p3d ax) ax)))))
756                      (let* ((t (list->kd-tree ax))
757                             (n (kd-tree-size t)))
758
759                        (list t boundary))))
760                ))
761            ))
762          )
763
764
765        (define (ClusteredRandomPointProcess cluster-pts n mean-distance x-seed y-seed boundary)
766
767          (match-let (((top left bottom right)
768                       (cases layer-boundary boundary
769                                   (Bounds (top left bottom right) 
770                                           (list top left bottom right))
771                                   (BoundsXZ (top left bottom right n k x y d d2) 
772                                             (list top left bottom right))
773                                   (BoundsYZ (top left bottom right n k x y d d2) 
774                                             (list top left bottom right))
775                                   )))
776
777
778          (let* (
779                 (x-extent   (- right left))
780                 (y-extent   (- top bottom))
781                 (z-extent-function (boundary-z-extent-function boundary))
782                 )
783
784            (let recur ((pts '()) (x-seed x-seed) (y-seed y-seed))
785
786              (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
787                    (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
788                    (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
789               
790                (let ((point-rejection1 (point2d-rejection top left bottom right)))
791                 
792                  (let inner-recur ((j 0) (ax pts))
793                    (if (< j n)
794                        (let ((x (f64vector-ref x-points j))
795                              (y (f64vector-ref y-points j))
796                              (z (f64vector-ref z-points j)))
797                          (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
798                            (let ((p3d 
799                                   (and (point-rejection1 p)
800                                        (compute-point3d p z z-extent-function))))
801                              (let ((pp (cluster-point-rejection p3d cluster-pts mean-distance)))
802                                (inner-recur (+ 1 j)  (if pp (cons pp ax) ax))))))
803
804                        (if (< (length ax) n)
805                            (recur ax (+ 1 x-seed) (+ 1 y-seed))
806                            (let* ((t (list->kd-tree (take ax n)))
807                                   (n (kd-tree-size t)))
808                             
809                              (list t boundary))))
810                    ))
811                ))
812            ))
813          )
814
815
816
817       
818        (define comment-pat (string->irregex "^#.*"))
819
820
821        (define (load-points-from-file filename)
822
823          (let ((in (open-input-file filename)))
824           
825            (if (not in) (error 'load-points-from-file "file not found" filename))
826
827            (let* ((lines
828                    (filter (lambda (line) (not (irregex-match comment-pat line)))
829                            (read-lines in)))
830
831                   (point-data
832                    (filter-map
833                     (lambda (line) 
834                       (let ((lst (map string->number (string-split line " \t"))))
835                         (and (not (null? lst)) (apply make-point lst))))
836                     lines))
837
838                   (point-tree (list->kd-tree point-data))
839                   )
840             
841              (list point-tree)
842             
843              ))
844          )
845
846
847        (define (segment-projection label source-tree target-sections zone my-comm myrank size)
848
849          (MPI:barrier my-comm)
850         
851          (let ((my-results
852                 (genpoint-projection label my-comm myrank size target-sections source-tree zone 0 0)))
853
854            (MPI:barrier my-comm)
855
856            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
857              (lambda (out-sources)
858                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
859                  (lambda (out-targets)
860                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
861                      (lambda (out-distances)
862                        (call-with-output-file (sprintf "~Asegments~A.dat"  label (if (> size 1) myrank ""))
863                          (lambda (out-segments)
864                            (for-each
865                             (lambda (my-data)
866                               (let* ((my-entry-len 5)
867                                      (my-data-len (/ (f64vector-length my-data) my-entry-len)))
868                                 (d "rank ~A: length my-data = ~A~%" myrank my-data-len)
869                                 (let recur ((m 0))
870                                   (if (< m my-data-len)
871                                       (let* (
872                                              (my-entry-offset (* m my-entry-len))
873                                              (source (inexact->exact (f64vector-ref my-data my-entry-offset)))
874                                              (target (inexact->exact (f64vector-ref my-data (+ 1 my-entry-offset))))
875                                              (distance (f64vector-ref my-data (+ 2 my-entry-offset)))
876                                              (segment (f64vector-ref my-data (+ 3 my-entry-offset)))
877                                              (section (f64vector-ref my-data (+ 4 my-entry-offset)))
878                                              )
879                                         (fprintf out-sources   "~A~%" source)
880                                         (fprintf out-targets   "~A~%" target)
881                                         (fprintf out-distances "~A~%" distance)
882                                         (fprintf out-segments  "~A ~A~%" segment section)
883                                         (recur (+ 1 m)))))
884                                 ))
885                             my-results))
886                          ))
887                      ))
888                  ))
889              ))
890          )
891       
892
893        (define (projection label source-tree target zone my-comm myrank size) 
894
895          (MPI:barrier my-comm)
896         
897          (let ((my-results (point-projection label my-comm myrank size target source-tree zone 0 (lambda (x nn) nn))))
898
899          (MPI:barrier my-comm)
900           
901            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
902              (lambda (out-sources)
903                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
904                  (lambda (out-targets)
905                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
906                      (lambda (out-distances)
907                        (for-each
908                         (lambda (my-data)
909                           (let* ((my-entry-len 3)
910                                  (my-data-len (/ (f64vector-length my-data) my-entry-len)))
911                             (d "~A: rank ~A: length my-data = ~A~%" label myrank my-data-len)
912                             (let recur ((m 0))
913                               (if (< m my-data-len)
914                                   (let ((source (inexact->exact (f64vector-ref my-data (* m my-entry-len))))
915                                         (target (inexact->exact (f64vector-ref my-data (+ 1 (* m my-entry-len)))))
916                                         (distance (f64vector-ref my-data (+ 2 (* m my-entry-len)))))
917                                     (fprintf out-sources "~A~%" source)
918                                     (fprintf out-targets "~A~%" target)
919                                     (fprintf out-distances "~A~%" distance)
920                                     (recur (+ 1 m)))))
921                             ))
922                         my-results)
923                        ))
924                    ))
925                ))
926            ))
927       
928
929        (include "calc-parser.scm")
930
931       
932        (define (load-config-file filename)
933          (let ((in (open-input-file filename)))
934            (if (not in) (error 'load-config-file "file not found" filename))
935            (init-bindings)
936            (let* ((lines (reverse (filter (lambda (x) (not (string-null? x))) (read-lines in))))
937                   (properties (filter-map
938                                (lambda (line) 
939                                  (and (not (string-prefix? "//" line))
940                                       (let ((lst (string-split line "=")))
941                                         (and (> (length lst) 1)
942                                              (let ((key (string->symbol (string-trim-both (car lst) #\space)))
943                                                    (val (calc-eval-string (cadr lst))))
944                                                (add-binding key val)
945                                                (cons key val))))))
946                                lines)))
947              properties
948              ))
949          )
950
951        (define (make-output-fname dirname sysname suffix . rest) 
952          (let-optionals 
953           rest ((x #t))
954           (and x
955                (let ((dirname (if (string? x) x dirname)))
956                  (let ((fname (sprintf "~A~A" sysname suffix)))
957                    (or (and dirname (make-pathname dirname fname)) fname)))
958                )))
959
960       
961)
962
963       
Note: See TracBrowser for help on using the repository browser.