source: project/release/4/npccl/trunk/npccl-utils.scm @ 30536

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

initial import of npccl, the Neural Parametric Curve Connectivity Language

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