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

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

picnic: additional imports in utils module

File size: 46.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 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 digraph graph-dfs)
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 concatenate)
41                (only srfi-4 
42                      s32vector s32vector-length s32vector-ref s32vector-set! make-s32vector
43                      f64vector f64vector? f64vector-ref f64vector-set! f64vector-length f64vector->list list->f64vector make-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        (define-record-type genpoint (make-genpoint coords parent-index parent-distance segment)
111          genpoint? 
112          (coords genpoint-coords)
113          (parent-index genpoint-parent-index)
114          (parent-distance genpoint-parent-distance)
115          (segment genpoint-segment)
116          )
117
118        (define-record-type swcpoint (make-swcpoint id type coords radius pre)
119          swcpoint? 
120          (id swcpoint-id)
121          (type swcpoint-type)
122          (coords swcpoint-coords)
123          (radius swcpoint-radius)
124          (pre swcpoint-pre)
125          )
126
127       
128        (define-record-type cell (make-cell ty index origin sections)
129          cell? 
130          (ty cell-type)
131          (index cell-index)
132          (origin cell-origin)
133          (sections cell-sections)
134          )
135       
136       
137        (define (cell-section-ref name cell)
138          (alist-ref name (cell-sections cell)))
139
140       
141        (define (write-pointset name pts)
142            (call-with-output-file (sprintf "~A.pointset.dat" name)
143              (lambda (out)
144                (for-each (match-lambda
145                           ((gid p)
146                            (fprintf out "~A ~A ~A~%" 
147                                     (coord 0 p)
148                                     (coord 1 p)
149                                     (coord 2 p))))
150                          pts))
151              ))
152
153       
154        (define (write-layout name pts #!optional rank)
155            (call-with-output-file (if rank 
156                                       (sprintf "~A.~A.layout.dat" name rank)
157                                       (sprintf "~A.layout.dat" name))
158              (lambda (out)
159                (for-each (match-lambda
160                           ((gid p)
161                            (fprintf out "~A ~A ~A ~A~%" 
162                                     gid
163                                     (coord 0 p)
164                                     (coord 1 p)
165                                     (coord 2 p))))
166                          pts))
167              ))
168
169        (define (write-sections forest-name section-name layout sections #!optional rank)
170          (call-with-output-file (if rank
171                                     (sprintf "~A.~A.~A.section.dat" forest-name section-name rank)
172                                     (sprintf "~A.~A.section.dat" forest-name section-name))
173            (lambda (out)
174              (for-each
175               (match-lambda*
176                (((gid p) section)
177                 (fprintf out "~A " gid)
178                 (for-each
179                  (lambda (neurites)
180                    (for-each
181                     (lambda (gd)
182                       (let ((p (genpoint-coords gd)))
183                         (fprintf out "~A ~A ~A "
184                                  (coord 0 p)
185                                  (coord 1 p)
186                                  (coord 2 p))))
187                     (cdr neurites)))
188                  (cdr section))
189                 (fprintf out "~%")
190                 ))
191               layout 
192               sections))))
193
194
195        (define (cells-sections->kd-tree cells section-name #!key
196                                         (make-value
197                                          (lambda (i v) 
198                                            (list (genpoint-parent-index v)
199                                                  (genpoint-parent-distance v))))
200                                         (make-point
201                                          (lambda (v) (genpoint-coords v))))
202          (let ((t 
203                 (let recur ((cells cells) (points '()))
204                   (if (null? cells)
205                       (list->kd-tree
206                        points
207                        make-value: make-value
208                        make-point: make-point)
209                       (let ((cell (car cells)))
210                         (recur (cdr cells) 
211                                (let inner ((sections (append (cell-section-ref section-name cell)))
212                                            (points points))
213                                  (if (null? sections) points
214                                      (inner
215                                       (cdr sections)
216                                       (append (cdr (car sections)) points))
217                                      ))
218                                ))
219                       ))
220                 ))
221            t))
222
223        (define (sections->kd-tree cells #!key
224                                   (make-value
225                                    (lambda (i v) 
226                                      (list (genpoint-parent-index v)
227                                            (genpoint-parent-distance v))))
228                                   (make-point
229                                    (lambda (v) (genpoint-coords v))))
230          (let ((t 
231                 (let recur ((cells cells) (points '()))
232                   (if (null? cells)
233                       (list->kd-tree
234                        points
235                        make-value: make-value
236                        make-point: make-point)
237                       (let ((sections (car cells)))
238                         (let inner ((sections sections) (points points))
239                           (if (null? sections)
240                               (recur (cdr cells) points)
241                               (let ((section (car sections)))
242                                 (inner (cdr sections) 
243                                        (append (cdr (cadr section)) points))
244                                 ))
245                           ))
246                       ))
247                 ))
248            t))
249
250
251        (define (cells-origins->kd-tree cells)
252          (let ((t 
253                 (let recur ((cells cells) (points '()))
254                   (if (null? cells)
255                       (list->kd-tree
256                        points
257                        make-value: (lambda (i v) 
258                                      (list (genpoint-parent-index v)
259                                            (genpoint-parent-distance v)))
260                        make-point: (lambda (v) (genpoint-coords v))
261                        )
262                       (let ((cell (car cells)))
263                         (recur (cdr cells) 
264                                (cons (make-genpoint (cell-origin cell) (cell-index cell) 0. 0)
265                                      points)))
266                       ))
267                 ))
268            t))
269
270        (define (pointCoord axis p)
271          (coord (inexact->exact axis) p))
272
273        ;; Samples a parametric curve at regular intervals in the range xmin..xmax inclusive.
274        (define (sample-uniform)
275          (lambda (n)
276            (let* (
277                   (xmin  0.0)
278                   (xmax  1.0)
279                   (delta (- xmax xmin))
280                   (dx    (if (zero? delta) 0
281                              (if (< n 2)
282                                (error 'sample-uniform "number of iterations must be >= 2")
283                                (/ (- xmax xmin) (- n 1)))))
284                   )
285              (list-tabulate n (lambda (i) (+ xmin (* dx i))))
286              )))
287
288
289        ;; Samples a parametric curve according to a polynomial
290        ;; c0 + c1 x + c2 x^2 in the range xmin..xmax inclusive.
291        (define (sample-polynomial c0 c1 c2)
292          (lambda (n)
293            (let* (
294                   (f (lambda (x) (+ c0 (* x c1) (* x x c2))))
295                   (dx (/ 1.0 n))
296                   )
297              (let recur ((i 0.0) (fi 0.0) (lst '()))
298                (let ((fi1 (+ fi (f i))))
299                  (if (<= fi1 1.0)
300                      (recur (+ i dx) fi1 (cons fi1 lst))
301                      (reverse lst))))
302              )))
303
304
305        (define (compose-curves c1 c2)
306          (let ((c (pcurve:compose-curve (list + + +) c1 c2)))
307;            (print "compose-curves: c1 = ") (pp (pcurve:iterate-curve c1 10))
308;            (print "compose-curves: c2 = ") (pp (pcurve:iterate-curve c2 10))
309;            (print "compose-curves: c = ") (pp (pcurve:iterate-curve c 10))
310            c))
311             
312
313        (define (make-simple-curve fx fy fz n) 
314          (let ((c (pcurve:simple-curve n 1 (list fx fy fz) 0.0 1.0)))
315            c
316            ))
317       
318
319        (define (make-harmonic axis amp period phase n)
320          (let* ((freq (/ (* 2 PI) (/ 1.0 period)))
321                 (pf   (lambda (t) (* amp (cos (+ (* freq t) phase)))))
322                 (zf   (lambda (t) 0.0))
323                 (c (pcurve:simple-curve 
324                    (inexact->exact n) 1
325                    (cond ((= axis 0)
326                           (list pf zf zf))
327                          ((= axis 1)
328                           (list zf pf zf))
329                          ((= axis 2)
330                           (list pf zf zf))
331                          (else
332                           (error 'make-harmonic "invalid axis" axis))
333                          )
334                    0.0 1.0)))
335            c
336            ))
337       
338
339        (define (make-line-segment p dx dy dz) 
340          (let ((c (pcurve:line-segment 3 (list dx dy dz))))
341            (pcurve:translate-curve (point->list p) c)
342            ))
343       
344       
345        ;; auxiliary function to create segment points
346        (define (make-segment si np sp xyz)
347          (list-tabulate 
348           np
349           (lambda (i) (let* ((pi (+ sp i))
350                              (p (make-point 
351                                  (f64vector-ref (car xyz) pi)
352                                  (f64vector-ref (cadr xyz) pi)
353                                  (f64vector-ref (caddr xyz) pi))))
354                         (list si p)
355                         ))
356           ))
357
358
359         
360
361        ;;
362        ;; Creates a process of the given segments and number of points per
363        ;; segment from the given curve.
364        ;;
365        (define (make-segmented-process c f ns np)
366          (let ((xyz ((pcurve:sample-curve* c) (f (* ns np)))))
367            (let ((np1 (floor (/ (f64vector-length (car xyz)) ns))))
368              (let recur ((si ns) (ax '()))
369                (if (positive? si)
370                    (let ((sp (* (- si 1) np1)))
371                      (recur (- si 1) (append (make-segment si np1 sp xyz) ax)))
372                    ax)
373                ))
374            ))
375
376        ;;
377        ;; Non-segmented process
378        ;;
379        (define (make-process c f np)
380          (let ((xyz ((pcurve:sample-curve* c) (f np))))
381            (list-tabulate 
382             (f64vector-length (car xyz))
383             (lambda (i) 
384               (make-point 
385                (f64vector-ref (car xyz) i)
386                (f64vector-ref (cadr xyz) i)
387                (f64vector-ref (caddr xyz) i)))
388             ))
389          )
390       
391       
392        ;;
393        ;; Creates a named section containing points from the given segmented processes.
394        ;;
395        (define (make-segmented-section gid p label ps)
396          `(,label . 
397                   ,(fold (lambda (i+proc ax)
398                            (let ((seci (car i+proc)) 
399                                  (proc (cadr i+proc)))
400                              (cons
401                               `(,seci . 
402                                       ,(map (lambda (sp)
403                                               (let ((segi (car sp))
404                                                     (dpoint (cadr sp)))
405                                                 (let ((soma-distance (sqrt (dist2 p dpoint))))
406                                                   (make-genpoint dpoint gid soma-distance segi))
407                                                 ))
408                                             proc))
409                               ax)
410                              ))
411                          '() ps)
412                   ))
413       
414        ;;
415        ;; Non-segmented named section
416        ;;
417        (define (make-section gid p label ps)
418          `(,label . 
419                   ,(fold (lambda (i+proc ax)
420                            (let* ((seci (car i+proc)) 
421                                   (proc (cadr i+proc))
422                                   (pts (map (lambda (dpoint)
423                                               (let ((soma-distance (sqrt (dist2 p dpoint))))
424                                                 (make-genpoint dpoint gid soma-distance #f)))
425                                             proc)))
426                              (d "make-section: label = ~A gid = ~A p = ~A pts = ~A~%" 
427                                 label gid p proc)
428                              (cons `(,seci . ,pts) ax)
429                              ))
430                          '() ps)
431                   ))
432       
433       
434       
435        (define (make-gen-kd-tree data #!key (threshold 0.0))
436         
437          (define (update-bb p x-min y-min z-min x-max y-max z-max)
438            (let ((x (coord 0 p)) (y (coord 1 p)) (z (coord 2 p)))
439              (if (< x (x-min)) (x-min x))
440              (if (< y (y-min)) (y-min y))
441              (if (< z (z-min)) (z-min z))
442              (if (> x (x-max)) (x-max x))
443              (if (> y (y-max)) (y-max y))
444              (if (> z (z-max)) (z-max z))
445              ))
446
447
448          (let* (
449                 (t (list->kd-tree
450                     (filter (lambda (x) (<= threshold (cdr x))) data)
451                     make-value: (lambda (i v) (cdr v))
452                     make-point: (lambda (v) (apply make-point (car v)))
453                     offset: 0
454                     ))
455                 (n (kd-tree-size t))
456                 (x-min (make-parameter +inf.0))
457                 (y-min (make-parameter +inf.0))
458                 (z-min (make-parameter +inf.0))
459                 (x-max (make-parameter -inf.0))
460                 (y-max (make-parameter -inf.0))
461                 (z-max (make-parameter -inf.0))
462                 (bb (begin (kd-tree-for-each
463                             (lambda (p) (update-bb p x-min y-min z-min
464                                                    x-max y-max z-max))
465                             t)
466                            (list (x-min) (y-min) (z-min) (x-max) (y-max) (z-max))))
467                 )
468
469            (cons bb t)
470
471            ))
472
473
474
475
476
477        (define (genpoint-projection prefix my-comm myrank size cells fibers zone cell-start fiber-start)
478
479          (d "rank ~A: prefix = ~A zone = ~A~%" myrank prefix zone)
480
481          (fold (lambda (cell ax)
482
483                  (let* ((i (+ cell-start (car cell)))
484                         (root (modulo i size))
485                         (sections (cadr cell)))
486                   
487                    (fold 
488                     
489                     (lambda (sec ax)
490                       
491                       (let ((seci (car sec))
492                             (gxs  (cdr sec)))
493                         
494                         (let ((query-data
495                                (fold 
496                                 (lambda (gd ax)
497                                   (d "rank ~A: querying point ~A (~A)~%" 
498                                      myrank (genpoint-coords gd) 
499                                      (genpoint-parent-index gd))
500                                   (fold
501                                    (lambda (x ax) 
502                                      (let (
503                                            (source (car x))
504                                            (target i)
505                                            (distance (cadr x))
506                                            (segi (genpoint-segment gd))
507                                            )
508                                        (if (not segi)
509                                            (error 'genpoint-projection "point does not have segment index" gd))
510                                        (append (list source target distance segi seci) ax)
511                                        ))
512                                    ax
513                                    (delete-duplicates
514                                     (map (lambda (x) 
515                                            (d "rank ~A: query result = ~A (~A) (~A) ~%" 
516                                               myrank (kdnn-point x) (kdnn-distance x) (kdnn-parent-index x))
517                                            (list (+ fiber-start (kdnn-parent-index x))
518                                                  (+ (kdnn-distance x) (genpoint-parent-distance gd)  (kdnn-parent-distance x))
519                                                  ))
520                                          (kd-tree-near-neighbors* fibers zone (genpoint-coords gd)))
521                                     (lambda (u v)
522                                       (= (car u) (car v)))
523                                     )
524                                    ))
525                                 '() gxs)))
526                           (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
527
528                                  (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
529                             
530                             (append res1 ax))
531                           
532                           ))
533                       )
534                     ax sections)
535                    ))
536                '() cells ))
537
538
539       
540
541        (define (point-projection prefix my-comm myrank size pts fibers zone point-start nn-filter)
542          (fold (lambda (px ax)
543
544                  (d "~A: rank ~A: px = ~A zone=~A~%"  prefix myrank px zone)
545
546                  (let* ((i (+ point-start (car px)))
547                         (root (modulo i size))
548                         (dd (d "~A: rank ~A: querying point ~A (root ~A)~%" prefix myrank px root))
549                         (query-data
550                          (let ((pd (cadr px))) 
551                            (fold
552                             (lambda (x ax) 
553                               (let ((source (car x))
554                                     (target i)
555                                     (distance (cadr x)))
556                                 (if (and (> distance  0.) (not (= source target)))
557                                     (append (list source target distance) ax)
558                                     ax)
559                                 ))
560                             '()
561                             (delete-duplicates
562                              (map (lambda (x) 
563                                     (let ((res (list (car (cadar x)) 
564                                                      (+ (caddr x) (cadr (cadar x))))))
565                                       (d "~A: x = ~A res = ~A~%" prefix x res)
566                                       res))
567                                   (nn-filter pd (kd-tree-near-neighbors* fibers zone pd))
568                                   )
569                              (lambda (u v) (= (car u) (car v)))
570                              )
571                             ))
572                          ))
573
574
575                    (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
576                           (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
577                      (append res1 ax))
578                    ))
579
580                '() pts))
581       
582
583
584
585
586        (define-datatype layer-boundary layer-boundary?
587          (Bounds (top number?) (left number?) (bottom number?) (right number?))
588          (BoundsXZ (top number?) (left number?) (bottom number?) (right number?)
589                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
590          (BoundsYZ (top number?) (left number?) (bottom number?) (right number?)
591                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
592          )
593
594
595
596        (define (boundary-z-extent-function boundary)
597          (cases layer-boundary boundary
598                 (Bounds (top left bottom right) 
599                         (lambda (x y) 0.))
600                 (BoundsXZ (top left bottom right n k x y d d2) 
601                           (lambda (xp yp) 
602                             (let-values (((y0tab y1tab y2tab res)
603                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
604                               (f64vector-ref y0tab 0))))
605                 (BoundsYZ (top left bottom right n k x y d d2) 
606                           (lambda (xp yp) 
607                             (let-values (((y0tab y1tab y2tab res)
608                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
609                               (f64vector-ref y0tab 0))))
610                 ))
611
612
613        (define (point2d-rejection top left bottom right)
614            (lambda (p)
615              (let ((x (coord 0 p)) (y (coord 1 p)))
616                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
617            )
618
619
620        (define (compute-point3d p zu z-extent-function)
621          (let* ((x (coord 0 p))
622                 (y (coord 1 p))
623                 (z-extent (z-extent-function x y)))
624            (if (zero? zu)
625                (make-point x y zu)
626                (make-point x y (fp* zu z-extent)))
627            ))
628
629
630        (define (cluster-point-rejection p cluster-pts mean-distance)
631          (let ((D (* 4 mean-distance mean-distance))
632                (nn (kd-tree-nearest-neighbor cluster-pts p)))
633            (and (< (dist2 p nn) D) p)))
634
635
636
637        (define (XZAxis n k x-points z-points boundary)
638         
639          (if (not (>= n 3)) 
640              (error 'generate-boundary "XZAxis requires at least 3 interpolation points (n >= 3)"))
641         
642          (cases layer-boundary boundary
643                 (Bounds (top left bottom right) 
644
645                         (let-values (((d d2 constr errc diagn)
646                                       (bvsp-spline:compute n k x-points z-points)))
647                           
648                           (if (not (zero? errc)) 
649                               (error 'generate-boundary "error in constructing spline from boundary points" errc))
650                           
651                           (BoundsXZ top left bottom right n k x-points z-points d d2)))
652                 
653                 (else (error 'generate-boundary "boundary argument to XZAxis is already a pseudo-3D boundary")))
654          )
655
656
657        (define (Grid x-spacing y-spacing z-spacing boundary)
658
659          (match-let (((top left bottom right)
660                       (cases layer-boundary boundary
661                                   (Bounds (top left bottom right) 
662                                           (list top left bottom right))
663                                   (BoundsXZ (top left bottom right n k x y d d2) 
664                                             (list top left bottom right))
665                                   (BoundsYZ (top left bottom right n k x y d d2) 
666                                             (list top left bottom right))
667                                   )))
668
669          (let* (
670                 (x-extent   (- right left))
671                 (y-extent   (- top bottom))
672                 (z-extent-function
673                  (boundary-z-extent-function boundary))
674                 (compute-grid-points3d
675                  (lambda (p z-spacing z-extent-function)
676                    (let* ((x (coord 0 p))
677                           (y (coord 1 p))
678                           (z-extent (z-extent-function x y)))
679                      (let recur ((points '()) (z (/ z-spacing 2.)))
680                        (if (> z z-extent)
681                            points
682                            (recur (cons (make-point x y z) points) (+ z z-spacing)))
683                        ))
684                    ))
685                 )
686           
687            (d "Grid: x-spacing = ~A~%" x-spacing)
688            (d "Grid: y-spacing = ~A~%" y-spacing)
689            (d "Grid: z-spacing = ~A~%" z-spacing)
690           
691            (d "Grid: x-extent = ~A~%" x-extent)
692            (d "Grid: y-extent = ~A~%" y-extent)
693           
694            (let ((x-points (let recur ((points '()) (x (/ x-spacing 2.)))
695                              (if (> x x-extent)
696                                  (list->f64vector (reverse points))
697                                  (recur (cons x points) (+ x x-spacing)))))
698                  (y-points (let recur ((points '()) (y (/ y-spacing 2.)))
699                              (if (> y y-extent)
700                                  (list->f64vector (reverse points))
701                                  (recur (cons y points) (+ y y-spacing)))))
702                  )
703             
704              (let ((nx (f64vector-length x-points))
705                    (ny (f64vector-length y-points))
706                    )
707               
708                (let recur ((i 0) (j 0) (ax '()))
709                  (if (< i nx)
710                      (let ((x (f64vector-ref x-points i)))
711                        (if (< j ny)
712                            (let* ((y   (f64vector-ref y-points j))
713                                   (p   (make-point x y))
714                                   (p3ds (if (zero? z-spacing)
715                                             (list (make-point x y 0.0))
716                                             (compute-grid-points3d p z-spacing z-extent-function)))
717                                   )
718                              (recur i (+ 1 j) (if p3ds (append p3ds ax) ax)))
719                            (recur (+ 1 i) 0 ax)))
720                      (let* ((t (list->kd-tree ax))
721                             (n (kd-tree-size t)))
722                        (list t boundary)
723                        ))
724                  )))
725            ))
726          )
727
728        (define (UniformRandomPointProcess n x-seed y-seed boundary)
729
730          (match-let (((top left bottom right)
731                       (cases layer-boundary boundary
732                                   (Bounds (top left bottom right) 
733                                           (list top left bottom right))
734                                   (BoundsXZ (top left bottom right n k x y d d2) 
735                                             (list top left bottom right))
736                                   (BoundsYZ (top left bottom right n k x y d d2) 
737                                             (list top left bottom right))
738                                   )))
739
740          (let* (
741                 (x-extent   (- right left))
742                 (y-extent   (- top bottom))
743                 (z-extent-function (boundary-z-extent-function boundary))
744                 )
745
746            (let ((x-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init x-seed)))
747                  (y-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init y-seed)))
748                  (z-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init (current-milliseconds)))))
749             
750              (let ((point-rejection1 (point2d-rejection top left bottom right)))
751               
752                (let recur ((i 0) (ax '()))
753                  (if (< i n)
754                      (let ((x (f64vector-ref x-points i))
755                            (y (f64vector-ref y-points i))
756                            (z (f64vector-ref z-points i)))
757                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
758                          (let ((p3d 
759                                 (and (point-rejection1 p)
760                                      (compute-point3d p z z-extent-function))))
761                            (recur (+ 1 i) (if p3d (cons p3d ax) ax)))))
762                      (let* ((t (list->kd-tree ax))
763                             (n (kd-tree-size t)))
764
765                        (list t boundary))))
766                ))
767            ))
768          )
769
770
771        (define (ClusteredRandomPointProcess cluster-pts n mean-distance x-seed y-seed boundary)
772
773          (match-let (((top left bottom right)
774                       (cases layer-boundary boundary
775                                   (Bounds (top left bottom right) 
776                                           (list top left bottom right))
777                                   (BoundsXZ (top left bottom right n k x y d d2) 
778                                             (list top left bottom right))
779                                   (BoundsYZ (top left bottom right n k x y d d2) 
780                                             (list top left bottom right))
781                                   )))
782
783
784          (let* (
785                 (x-extent   (- right left))
786                 (y-extent   (- top bottom))
787                 (z-extent-function (boundary-z-extent-function boundary))
788                 )
789
790            (let recur ((pts '()) (x-seed x-seed) (y-seed y-seed))
791
792              (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
793                    (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
794                    (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
795               
796                (let ((point-rejection1 (point2d-rejection top left bottom right)))
797                 
798                  (let inner-recur ((j 0) (ax pts))
799                    (if (< j n)
800                        (let ((x (f64vector-ref x-points j))
801                              (y (f64vector-ref y-points j))
802                              (z (f64vector-ref z-points j)))
803                          (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
804                            (let ((p3d 
805                                   (and (point-rejection1 p)
806                                        (compute-point3d p z z-extent-function))))
807                              (let ((pp (cluster-point-rejection p3d cluster-pts mean-distance)))
808                                (inner-recur (+ 1 j)  (if pp (cons pp ax) ax))))))
809
810                        (if (< (length ax) n)
811                            (recur ax (+ 1 x-seed) (+ 1 y-seed))
812                            (let* ((t (list->kd-tree (take ax n)))
813                                   (n (kd-tree-size t)))
814                             
815                              (list t boundary))))
816                    ))
817                ))
818            ))
819          )
820
821
822
823       
824        (define comment-pat (string->irregex "^#.*"))
825
826
827        (define (load-points-from-file filename)
828
829          (let ((in (open-input-file filename)))
830           
831            (if (not in) (error 'load-points-from-file "file not found" filename))
832
833            (let* ((lines
834                    (filter (lambda (line) (not (irregex-match comment-pat line)))
835                            (read-lines in)))
836
837                   (point-data
838                    (filter-map
839                     (lambda (line) 
840                       (let ((lst (map string->number (string-split line " \t"))))
841                         (and (not (null? lst)) (apply make-point lst))))
842                     lines))
843
844                   (point-tree (list->kd-tree point-data))
845                   )
846             
847              (list point-tree)
848             
849              ))
850          )
851
852
853        (define (make-tree-graph lst label)
854         
855          (let* (
856                 (g              (make-digraph label #f))
857                 (node-info      (g 'node-info))
858                 (node-info-set! (g 'node-info-set!))
859                 (add-node!      (g 'add-node!))
860                 (add-edge!      (g 'add-edge!))
861                 )
862
863            ;; insert nodes
864            (let recur ((lst lst))
865
866              (if (not (null? lst))
867
868                  (let ((point (car lst)))
869
870                    (let ((node-id (swcpoint-id point)))
871
872                      (add-node! node-id point)
873
874                      (recur (cdr lst))))))
875
876            ;; insert edges
877            (let recur ((lst lst))
878             
879              (if (not (null? lst))
880                 
881                  (let ((point (car lst)))
882                   
883                    (let ((node-id (swcpoint-id point))
884                          (pre-id  (swcpoint-pre point)))
885                     
886                      (let* ((pre-point   (node-info pre-id))
887                             (pre-coords  (swcpoint-coords pre-point))
888                             (node-coords (swcpoint-coords point))
889                             (distance    (sqrt (dist2 node-coords pre-coords))))
890                       
891                        (add-edge! (list node-id pre-id distance))
892                       
893                        (recur (cdr lst))
894                        ))
895                    ))
896              )
897            g 
898            ))
899
900
901        (define (tree-graph-distances+segments g nseg)
902
903          (define n        ((g 'capacity)))
904          (define distv    (make-f64vector (+ 1 n) -1.0))
905          (define rdistv   (make-f64vector (+ 1 n) -1.0))
906          (define segv     (make-s32vector (+ 1 n) -1))
907
908
909          ;; determine distances from root
910          (define (traverse-dist es)
911            (if (null? es) distv
912                (match-let (((i j dist) (car es)))
913                  (if (>= (f64vector-ref distv j) 0.0)
914                      (traverse-dist (cdr es))
915                      (let ((idist (f64vector-ref distv i)))
916                        (f64vector-set! distv j (+ idist dist))
917                        (let ((distv1 (traverse-dist ((g 'out-edges) j))))
918                          (traverse-dist es)))
919                      ))
920                ))
921         
922         
923          ;; determine distances from end (reverse distance)
924          (define (traverse-rdist es)
925            (if (null? es) rdistv
926                (match-let (((i j dist) (car es)))
927                  (if (>= (f64vector-ref rdistv i) 0.0)
928                      (traverse-rdist (cdr es))
929                      (let ((jdist (f64vector-ref distv j)))
930                        (f64vector-set! rdistv i (+ jdist dist))
931                        (let ((rdistv1 (traverse-rdist ((g 'in-edges) i))))
932                          (traverse-rdist es)))
933                      ))
934                ))
935
936
937          (define (compute-segv distv rdistv)
938            (let recur ((n n))
939              (if (>= n 1)
940                  (let* ((dist  (f64vector-ref distv n))
941                         (rdist (f64vector-ref rdistv n))
942                         (len   (+ dist rdist))
943                         (delta (round (/ len nseg)))
944                         (seg   (round (/ dist delta))))
945                    (s32vector-set! segv n (exact->inexact seg))
946                    (recur (- n 1))
947                    ))
948              ))
949         
950          (let ((root 1) 
951                (in-edges (g 'in-edges)) 
952                (terminals ((g 'terminals))))
953            (f64vector-set! distv root 0.0)
954            (for-each (lambda (x) (f64vector-set! distv x 0.0)) terminals)
955            (s32vector-set! segv root 0)
956            (traverse-dist ((g 'out-edges) root))
957            (traverse-rdist (concatenate (map (lambda (x) (in-edges x)) terminals)))
958            (compute-segv distv rdistv)
959            (list distv segv)
960          ))
961
962
963        (define (tree-graph->genpoints g gdistv gsegv type cell-index)
964         
965          (let ((node-info (g 'node-info))
966                (out-edges (g 'out-edges)))
967
968            (let recur ((n 1) (lst '()))
969
970              (let* (
971                     (point (node-info n))
972                     (point-type (swcpoint-type point))
973                     (point-pre (swcpoint-pre point))
974                     (proceed? (or (= point-type type)
975                                   (case (swcpoint-type point)
976                                     ((0 1 5 6) #t)
977                                     (else #f))))
978                     )
979                 
980                (if proceed?
981
982                    (let (
983                          (point1 (make-genpoint 
984                                   (swcpoint-coords point)
985                                   cell-index 
986                                   (f64vector-ref gdistv n)
987                                   (s32vector-ref gsegv n)))
988                          )
989
990                      (fold (lambda (x ax) (recur x ax)) (cons point1 lst) (out-edges n)))
991
992                    lst)
993
994                ))
995            ))
996
997 
998        (define (load-swc filename label type nseg cell-index)
999         
1000          (let ((in (open-input-file filename)))
1001           
1002            (if (not in) (error 'load-swc "file not found" filename))
1003           
1004            (let* (
1005                   (lines
1006                    (filter (lambda (line) (not (irregex-match comment-pat line)))
1007                            (read-lines in)))
1008
1009                   (swc-data
1010                    (filter-map
1011                     (lambda (line) 
1012                       (let ((lst (map string->number (string-split line " \t"))))
1013                         (and (not (null? lst)) 
1014                              (match-let (((id my-type x y z radius pre) lst))
1015                                         (make-swcpoint id type (make-point x y z)
1016                                                        radius pre)))
1017                         ))
1018                     lines))
1019
1020                   (swc-graph (make-tree-graph swc-data label))
1021
1022                   (dist+segs  (tree-graph-distances+segments swc-graph nseg))
1023                   (point-data (tree-graph->genpoints 
1024                                swc-graph (car dist+segs) (cadr dist+segs) 
1025                                type cell-index))
1026
1027                   )
1028
1029              (list point-data))
1030          ))
1031
1032
1033
1034
1035        (define (segment-projection label source-tree target-sections zone my-comm myrank size)
1036
1037          (MPI:barrier my-comm)
1038         
1039          (let ((my-results
1040                 (genpoint-projection label my-comm myrank size target-sections source-tree zone 0 0)))
1041
1042            (MPI:barrier my-comm)
1043
1044            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
1045              (lambda (out-sources)
1046                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
1047                  (lambda (out-targets)
1048                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
1049                      (lambda (out-distances)
1050                        (call-with-output-file (sprintf "~Asegments~A.dat"  label (if (> size 1) myrank ""))
1051                          (lambda (out-segments)
1052                            (for-each
1053                             (lambda (my-data)
1054                               (let* ((my-entry-len 5)
1055                                      (my-data-len (/ (f64vector-length my-data) my-entry-len)))
1056                                 (d "rank ~A: length my-data = ~A~%" myrank my-data-len)
1057                                 (let recur ((m 0))
1058                                   (if (< m my-data-len)
1059                                       (let* (
1060                                              (my-entry-offset (* m my-entry-len))
1061                                              (source (inexact->exact (f64vector-ref my-data my-entry-offset)))
1062                                              (target (inexact->exact (f64vector-ref my-data (+ 1 my-entry-offset))))
1063                                              (distance (f64vector-ref my-data (+ 2 my-entry-offset)))
1064                                              (segment (f64vector-ref my-data (+ 3 my-entry-offset)))
1065                                              (section (f64vector-ref my-data (+ 4 my-entry-offset)))
1066                                              )
1067                                         (fprintf out-sources   "~A~%" source)
1068                                         (fprintf out-targets   "~A~%" target)
1069                                         (fprintf out-distances "~A~%" distance)
1070                                         (fprintf out-segments  "~A ~A~%" segment section)
1071                                         (recur (+ 1 m)))))
1072                                 ))
1073                             my-results))
1074                          ))
1075                      ))
1076                  ))
1077              ))
1078          )
1079       
1080
1081        (define (projection label source-tree target zone my-comm myrank size) 
1082
1083          (MPI:barrier my-comm)
1084         
1085          (let ((my-results (point-projection label my-comm myrank size target source-tree zone 0 (lambda (x nn) nn))))
1086
1087          (MPI:barrier my-comm)
1088           
1089            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
1090              (lambda (out-sources)
1091                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
1092                  (lambda (out-targets)
1093                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
1094                      (lambda (out-distances)
1095                        (for-each
1096                         (lambda (my-data)
1097                           (let* ((my-entry-len 3)
1098                                  (my-data-len (/ (f64vector-length my-data) my-entry-len)))
1099                             (d "~A: rank ~A: length my-data = ~A~%" label myrank my-data-len)
1100                             (let recur ((m 0))
1101                               (if (< m my-data-len)
1102                                   (let ((source (inexact->exact (f64vector-ref my-data (* m my-entry-len))))
1103                                         (target (inexact->exact (f64vector-ref my-data (+ 1 (* m my-entry-len)))))
1104                                         (distance (f64vector-ref my-data (+ 2 (* m my-entry-len)))))
1105                                     (fprintf out-sources "~A~%" source)
1106                                     (fprintf out-targets "~A~%" target)
1107                                     (fprintf out-distances "~A~%" distance)
1108                                     (recur (+ 1 m)))))
1109                             ))
1110                         my-results)
1111                        ))
1112                    ))
1113                ))
1114            ))
1115       
1116
1117        (include "calc-parser.scm")
1118
1119       
1120        (define (load-config-file filename)
1121          (let ((in (open-input-file filename)))
1122            (if (not in) (error 'load-config-file "file not found" filename))
1123            (init-bindings)
1124            (let* ((lines (reverse (filter (lambda (x) (not (string-null? x))) (read-lines in))))
1125                   (properties (filter-map
1126                                (lambda (line) 
1127                                  (and (not (string-prefix? "//" line))
1128                                       (let ((lst (string-split line "=")))
1129                                         (and (> (length lst) 1)
1130                                              (let ((key (string->symbol (string-trim-both (car lst) #\space)))
1131                                                    (val (calc-eval-string (cadr lst))))
1132                                                (add-binding key val)
1133                                                (cons key val))))))
1134                                lines)))
1135              properties
1136              ))
1137          )
1138
1139        (define (make-output-fname dirname sysname suffix . rest) 
1140          (let-optionals 
1141           rest ((x #t))
1142           (and x
1143                (let ((dirname (if (string? x) x dirname)))
1144                  (let ((fname (sprintf "~A~A" sysname suffix)))
1145                    (or (and dirname (make-pathname dirname fname)) fname)))
1146                )))
1147
1148       
1149)
1150
1151       
Note: See TracBrowser for help on using the repository browser.