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

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

npccl: adding bounds constructor and random uniform point generation to the expression language

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