File: matrix.lisp

package info (click to toggle)
cl-graph 20180131-1
  • links: PTS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 240 kB
  • sloc: lisp: 2,884; makefile: 16
file content (714 lines) | stat: -rw-r--r-- 27,163 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
;;; graph-matrix.lisp --- build and manipulate matrix graph representations

;; Copyright (C) Eric Schulte and Tom Dye 2013

;; Licensed under the Gnu Public License Version 3 or later

;;; Commentary

;; Functions for manipulating matrix graph representations.

;;; Code:
(defpackage #:graph/matrix
  (:use :common-lisp
        :alexandria
        :metabang-bind
        :named-readtables
        :curry-compose-reader-macros
        :graph
        :fl.function)
  ;; shadow functions defined in alexandria, fl.function, and graph
  (:shadow :copy :factorial :standard-deviation :variance :median :mean :degree)
  (:export
   :matrix
   :fast-matrix
   :matrix-ref
   :matrix-n-rows
   :matrix-n-cols
   :matrix-same-size-p
   :matrix-symmetric-p
   :matrix-entries-different-p
   :matrix-copy
   :matrix-transpose
   :make-universal-matrix
   :make-identity-matrix
   :make-zeros-matrix
   :to-adjacency-matrix
   :to-reachability-matrix
   :reachablep
   :reachable-from
   :to-strong-component-matrix
   :strong-component-of
   :to-distance-matrix
   :distance-from-to
   :reflexivep
   :irreflexivep
   :symmetricp
   :asymmetricp
   :transitivep
   :intransitivep
   :completep
   :relational-structure
   :infinite
   :infinitep))
(in-package :graph/matrix)
(in-readtable :curry-compose-reader-macros)

(defclass matrix ()
  ((self :initarg :self :accessor self :initform nil)))

(defclass fast-matrix (matrix) ())

(defgeneric infinite (matrix)
  (:documentation "Return the most-positive value for the element type
  of MATRIX."))

(defmethod infinite ((matrix matrix))
  most-positive-fixnum)

(defmethod infinite ((matrix fast-matrix))
  most-positive-single-float)

(defgeneric infinitep (value matrix)
  (:documentation "Non-nil if VALUE is the most-positive value that
  can be held in MATRIX."))

(defmethod infinitep (value (matrix matrix))
  (= value (infinite matrix)))

(defgeneric matrix-ref (matrix row col)
  (:documentation "Return the value at ROW and COL in MATRIX."))

(defmethod matrix-ref ((matrix matrix) row col)
  (aref (self matrix) row col))

(defmethod matrix-ref ((fm fast-matrix) row col)
  (fl.function::mref (self fm) row col))

(defgeneric (setf matrix-ref) (new matrix row col)
  (:documentation "Make matrix-ref setf-able."))

(defmethod (setf matrix-ref) (new (matrix matrix) row col)
  (setf (aref (self matrix) row col) new))

(defmethod (setf matrix-ref) (new (fm fast-matrix) row col)
  (setf (fl.function::mref (self fm) row col) new))

(defgeneric matrix-n-rows (matrix)
  (:documentation "Return the number of rows in MATRIX."))

(defmethod matrix-n-rows ((matrix matrix))
  (if (self matrix)
      (array-dimension (self matrix) 0)
      0))

(defmethod matrix-n-rows ((matrix fast-matrix))
  (if (self matrix)
      (fl.function::nrows (self matrix))
      0))

(defgeneric matrix-n-cols (matrix)
  (:documentation "Return the number of columns in MATRIX."))

(defmethod matrix-n-cols ((matrix matrix))
  (if (self matrix)
      (array-dimension (self matrix) 1)
      0))

(defmethod matrix-n-cols ((matrix fast-matrix))
  (if (self matrix)
      (fl.function::ncols (self matrix))
      0))

(defun matrix-same-size-p (m1 m2)
  "Return t if matrix M1 has the same number of rows and columns as
matrix M2, nil otherwise."
  (and (= (matrix-n-rows m1) (matrix-n-rows m2))
       (= (matrix-n-cols m1) (matrix-n-cols m2))))

(defun matrix-entries-different-p (m1 m2)
  "Returns nil if the entries in matrix M1 and matrix M2 do not differ
from one another.  Returns 1 if the sizes of matrix M1 and matrix M2
differ.  Otherwise, returns a list of lists containing discrepant
entries. "
  (let ((result))
    (if (matrix-same-size-p m1 m2)
        (let ((m (matrix-n-rows m1))
              (n (matrix-n-cols m1)))
          (loop :for i :from 0 :below m :do
             (loop :for j :from 0 :below n :do
                (unless (= (matrix-ref m1 i j)
                           (matrix-ref m2 i j))
                  (push (list i j) result))))
          (when result (reverse result)))
        (setf result 1))
    result))

(defun matrix-symmetric-p (matrix)
  "Return t if matrix MATRIX is symmetric, nil otherwise."
  (not (matrix-entries-different-p matrix (matrix-transpose matrix))))

(defgeneric matrix-copy (matrix)
  (:documentation "Return a copy of MATRIX."))

(defmethod matrix-copy ((matrix matrix))
  (let* ((m (matrix-n-rows matrix))
         (n (matrix-n-cols matrix))
         (result (make-zeros-matrix (make-instance 'matrix) m n)))
    (when (self matrix)
      (loop :for i :from 0 :below m :do
         (loop :for j :from 0 :below n :do
            (setf (matrix-ref result i j) (matrix-ref matrix i j)))))
    result))

(defmethod matrix-copy ((fm fast-matrix))
  (let ((result (make-instance 'fast-matrix)))
    (when (self fm)
      (setf (self result) (fl.function::copy (self fm))))
    result))

(defgeneric matrix-sum (m1 m2 &key boolean)
  (:documentation "Return the result of adding matrix M1 and matrix
  M2. M1 and M2 must be the same size. If BOOLEAN is non-nil, then use
  boolean arithmetic, where 1+1=1."))

(defmethod matrix-sum ((m1 matrix) (m2 matrix) &key boolean)
  (and (matrix-same-size-p m1 m2)
       (let* ((m (matrix-n-rows m1))
              (n (matrix-n-cols m1))
              (result (make-zeros-matrix (make-instance 'matrix) m n))
              (zero 0)
              (one 1))
         (declare (type fixnum zero))
         (declare (type fixnum one))
         (loop :for i :from 0 :below m :do
            (loop :for j :from 0 :below n :do
               (setf (matrix-ref result i j)
                     (if boolean
                         (if (> (+ (matrix-ref m1 i j)
                                   (matrix-ref m2 i j)) 0) one zero)
                         (+ (matrix-ref m1 i j)
                            (matrix-ref m2 i j))))))
         result)))

(defmethod matrix-sum ((m1 fast-matrix) (m2 fast-matrix) &key boolean)
  (when (matrix-same-size-p m1 m2)
    (let ((result (make-instance 'fast-matrix)))
      (setf (self result) (fl.function::m+ (self m1) (self m2)))
      (when boolean
        (let ((m (matrix-n-rows result))
              (n (matrix-n-cols result))
              (one 1.0s0))
          (declare (type single-float one))
          (loop :for i :from 0 :below m :do
             (loop :for j :from 0 :below n :do
                (if (> (matrix-ref result i j) 0)
                    (setf (matrix-ref result i j) one))))))
      result)))

(defgeneric matrix-difference (m1 m2)
  (:documentation "Return the result of subtracting M2 from M1. M1 and
  M2 must be the same size."))

(defmethod matrix-difference ((m1 matrix) (m2 matrix))
  (and (matrix-same-size-p m1 m2)
       (let ((result (matrix-copy m1))
             (m (matrix-n-rows m1))
             (n (matrix-n-cols m1)))
         (loop :for i :from 0 :below m :do
            (loop :for j :from 0 :below n :do
               (setf (matrix-ref result i j)
                     (- (matrix-ref result i j)
                        (matrix-ref m2 i j)))))
         result)))

(defgeneric matrix-elementwise-product (m1 m2 &key boolean)
  (:documentation "Return the result of multiplying the elements of
  matrix M1 and matrix M2. M1 and M2 must be the same size."))

(defmethod matrix-elementwise-product ((m1 matrix) (m2 matrix) &key boolean)
  (and (matrix-same-size-p m1 m2)
       (let ((result (matrix-copy m1))
             (m (matrix-n-rows m1))
             (n (matrix-n-cols m1)))
         (loop :for i :from 0 :below m :do
            (loop :for j :from 0 :below n :do
               (setf (matrix-ref result i j)
                     (if boolean
                         (if (* (matrix-ref result i j)
                                (matrix-ref m2 i j)) 1 0)
                         (* (matrix-ref result i j)
                            (matrix-ref m2 i j))))))
         result)))

(defgeneric matrix-product (m1 m2)
  (:documentation "Return the result of multiplying matrix M1 and
  matrix M2. The number of columns of M1 must equal the number of rows
  of M2."))

(defmethod matrix-product ((m1 matrix) (m2 matrix))
  (and (= (matrix-n-cols m1) (matrix-n-rows m2))
       (loop
          :with m = (matrix-n-rows m1)
          :with n = (matrix-n-cols m1)
          :with l = (matrix-n-cols m2)
          :with c = (make-zeros-matrix (make-instance 'matrix) m l)
          :for i :below m :do
          (loop :for k :below l :do
             (setf (matrix-ref c i k)
                   (loop :for j :below n
                      :sum (* (matrix-ref m1 i j)
                              (matrix-ref m2 j k)))))
          :finally (return c)))  )

(defmethod matrix-product ((m1 fast-matrix) (m2 fast-matrix))
  (and (= (matrix-n-cols m1) (matrix-n-rows m2))
       (let ((result (make-instance 'fast-matrix)))
         (setf (self result) (fl.function::m* (self m1) (self m2)))
         result)))

(defgeneric matrix-transpose (matrix)
  (:documentation "Return a new matrix that interchanges the rows and
  columns of MATRIX."))

(defmethod matrix-transpose ((matrix matrix))
  (let ((m (matrix-n-rows matrix))
        (n (matrix-n-cols matrix))
        (result (make-instance 'matrix)))
    (setf result (make-zeros-matrix result n m))
    (loop :for i :from 0 :below m :do
       (loop :for j :from 0 :below n :do
          (setf (matrix-ref result j i)
                (matrix-ref matrix i j))))
    result))

(defmethod matrix-transpose ((fm fast-matrix))
  (let ((result (make-instance 'fast-matrix)))
    (setf (self result) (fl.function::transpose (self fm)))
    result))

(defgeneric make-zeros-matrix (matrix rows cols)
  (:documentation "Return matrix MATRIX with ROWS rows and COLS
  columns of zeros."))

(defmethod make-zeros-matrix ((matrix matrix) rows cols)
  (setf (self matrix) (make-array (list rows cols)
                                  :element-type 'fixnum
                                  :initial-element 0))
  matrix)

(defmethod make-zeros-matrix ((fm fast-matrix) rows cols)
  (setf (self fm) (fl.function::zeros rows cols 'single-float))
  fm)

(defgeneric make-universal-matrix (matrix rows cols)
  (:documentation "Return a universal matrix with ROWS rows and COLS columns."))

(defmethod make-universal-matrix ((matrix matrix) rows cols)
  (setf (self matrix) (make-array (list rows cols)
                                  :element-type 'fixnum
                                  :initial-element 1))
  matrix)

(defmethod make-universal-matrix ((fm fast-matrix) rows cols)
  (setf (self fm) (fl.function::ones rows cols 'single-float))
  fm)

(defgeneric make-infinity-matrix (matrix rows cols)
  (:documentation "Return a matrix of ROWS rows and COLS cols with
  each entry set to infinity"))

(defmethod make-infinity-matrix ((matrix matrix) rows cols)
  (progn
    (setf (self matrix) (make-array (list rows cols)
                                    :element-type 'fixnum
                                    :initial-element (infinite matrix)))
    matrix))

(defmethod make-infinity-matrix ((fm fast-matrix) rows cols)
  (progn
    (setf (self fm) (fl.function::zeros rows cols 'single-float))
    (fl.function::fill! (self fm) (infinite fm))
    ;; (loop :for i :from 0 :below rows :do
    ;;    (loop :for j :from 0 :below cols :do
    ;;       (setf (matrix-ref fm i j) infinity)))
    fm))

(defgeneric make-identity-matrix (matrix order)
  (:documentation "Return an identity matrix of order ORDER."))

(defmethod make-identity-matrix ((matrix matrix) order)
  (setf matrix (make-zeros-matrix matrix order order))
  (loop :for i :from 0 :below order :do
     (setf (matrix-ref matrix i i) 1))
  matrix)

(defmethod make-identity-matrix ((fm fast-matrix) order)
  (setf (self fm) (fl.function::eye order order 'single-float))
  fm)

;; Adapted from
;; https://rosettacode.org/wiki/Matrix-exponentiation_operator#Common_Lisp

(defgeneric matrix-power (matrix exp)
  (:documentation "Raise MATRIX to the power EXP and return the result."))

(defmethod matrix-power ((matrix matrix) exp)
  (let ((m-rows (matrix-n-rows matrix)))
    (cond
      ((/= m-rows (matrix-n-cols matrix)) (error "Non-square matrix"))
      ((zerop exp) (make-identity-matrix matrix m-rows))
      ((= 1 exp) (matrix-copy matrix))
      ((zerop (mod exp 2)) (let ((me2 (matrix-power matrix (/ exp 2))))
                             (matrix-product me2 me2)))
      (t (let ((me2 (matrix-power matrix (/ (1- exp) 2))))
           (matrix-product matrix (matrix-product me2 me2)))))))

(defgeneric to-adjacency-matrix (graph matrix)
  (:documentation "Return the adjacency matrix of GRAPH."))

(defmethod to-adjacency-matrix ((graph graph) (matrix matrix))
  (let ((node-index-hash (make-hash-table))
        (counter -1))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (setf matrix (make-zeros-matrix matrix (+ counter 1) (+ counter 1)))
    (mapc (lambda-bind ((a b))
                  (setf (matrix-ref matrix
                                    (gethash a node-index-hash)
                                    (gethash b node-index-hash))
                        1)
                  (setf (matrix-ref matrix
                                    (gethash b node-index-hash)
                                    (gethash a node-index-hash))
                        1))
          (edges graph))
    matrix))

(defmethod to-adjacency-matrix ((graph digraph) (matrix matrix))
  (let ((node-index-hash (make-hash-table))
        (counter -1))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (setf matrix (make-zeros-matrix matrix (+ counter 1) (+ counter 1)))
    (mapc (lambda-bind ((a b))
                  (setf (matrix-ref matrix
                                    (gethash a node-index-hash)
                                    (gethash b node-index-hash))
                        1))
          (edges graph))
    matrix))

(defmethod to-adjacency-matrix ((graph graph) (matrix fast-matrix))
  (let ((node-index-hash (make-hash-table))
        (counter -1)
        (one 1.0s0))
    (declare (type single-float one))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (setf matrix (make-zeros-matrix matrix (+ counter 1) (+ counter 1)))
    (mapc (lambda-bind ((a b))
                  (setf (matrix-ref matrix
                                    (gethash a node-index-hash)
                                    (gethash b node-index-hash))
                        one)
                  (setf (matrix-ref matrix
                                    (gethash b node-index-hash)
                                    (gethash a node-index-hash))
                        one))
          (edges graph))
    matrix))

(defmethod to-adjacency-matrix ((graph digraph) (matrix fast-matrix))
  (let ((node-index-hash (make-hash-table))
        (counter -1)
        (one 1.0s0))
    (declare (type single-float one))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (setf matrix (make-zeros-matrix matrix (+ counter 1) (+ counter 1)))
    (mapc (lambda-bind ((a b))
                  (setf (matrix-ref matrix
                                    (gethash a node-index-hash)
                                    (gethash b node-index-hash))
                        one))
          (edges graph))
    matrix))

(defgeneric to-reachability-matrix (graph matrix &key limit)
  (:documentation "Return the reachability matrix of the graph GRAPH.
  With the optional argument LIMIT set to an integer in the range 2 to
  two less than the number of nodes in GRAPH, produces a limited
  reachability matrix with paths of length LIMIT or less."))

(defmethod to-reachability-matrix ((graph graph) (matrix matrix) &key limit)
  (let ((n (length (nodes graph))))
    (assert (or (not limit)
                (and (integerp limit) (> limit 1) (< limit (- n 1))))
            (limit)
            "~S must be an integer between 2 and ~S"
            limit (- n 2))
    (let* ((result (make-identity-matrix (make-instance 'matrix) n))
           (max-power (or limit (- n 1)))
           (adjacency (to-adjacency-matrix graph (make-instance 'matrix)))
           (adjacency-powers (matrix-copy adjacency)))
      (setf result (matrix-sum adjacency result :boolean t))
      (loop :for i :from 2 :to max-power :do
         (setf adjacency-powers (matrix-product adjacency-powers adjacency))
         (setf result (matrix-sum adjacency-powers result :boolean t)))
      result)))

(defmethod to-reachability-matrix ((graph graph) (matrix fast-matrix) &key limit)
  (let ((n (length (nodes graph))))
    (assert (or (not limit)
                (and (integerp limit) (> limit 1) (< limit (- n 1))))
            (limit)
            "~S must be an integer between 2 and ~S"
            limit (- n 2))
    (let* ((result (make-identity-matrix (make-instance 'fast-matrix) n))
           (max-power (or limit (- n 1)))
           (adjacency (to-adjacency-matrix graph (make-instance 'fast-matrix)))
           (adjacency-powers (matrix-copy adjacency)))
      (setf result (matrix-sum adjacency result :boolean t))
      (loop :for i :from 2 :to max-power :do
         (setf adjacency-powers (matrix-product adjacency-powers adjacency))
         (setf result (matrix-sum adjacency-powers result :boolean t)))
      result)))

(defgeneric reachablep (graph rd from to)
  (:documentation "Given a graph GRAPH and a reachability matrix RD,
  returns t if node TO is reachable from node FROM, nil otherwise."))

(defmethod reachablep ((graph graph) (rd matrix) from to)
  (let ((node-index-hash (make-hash-table))
        (counter -1))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (= 1 (matrix-ref rd (gethash from node-index-hash)
                     (gethash to node-index-hash)))))

(defgeneric reachable-from (graph rd from)
  (:documentation "Given a reachability matrix RD, return a list of
  the nodes in graph GRAPH reachable from node FROM."))

(defmethod reachable-from ((graph graph) (rd matrix) from)
  (let ((node-index-hash (make-hash-table))
        (counter -1)
        (result))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (maphash #'(lambda (k v)
                 (unless
                     (= 0 (matrix-ref rd (gethash from node-index-hash) v))
                   (push k result)))
             node-index-hash)
    (reverse result)))

(defgeneric to-strong-component-matrix (rd)
  (:documentation "Given a reachability matrix of a digraph, RD,
  return a matrix in which the strong component of GRAPH containing
  node_i is given by the entries of 1 in the ith row (or column)."))

(defmethod to-strong-component-matrix ((rd matrix))
  (matrix-elementwise-product rd (matrix-transpose rd)))

(defgeneric strong-component-of (node graph strong-components)
  (:documentation "Return a list of nodes from graph GRAPH in the
  strong component that contains node NODE, as given by the strong
  component matrix STRONG-COMPONENTS."))

(defmethod strong-component-of (node (graph graph) (strong-components matrix))
  (let ((node-index-hash (make-hash-table))
        (counter -1)
        (result))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (maphash #'(lambda (k v)
                 (unless (= 0 (matrix-ref
                                   strong-components
                                   (gethash node node-index-hash) v))
                   (push k result)))
             node-index-hash)
    (reverse result)))

(defgeneric to-distance-matrix (graph nd)
  (:documentation "Return the distance matrix ND of graph GRAPH."))

(defmethod to-distance-matrix ((graph graph) (nd matrix))
  (let* ((a (to-adjacency-matrix graph (make-instance 'matrix)))
         (a-power (to-adjacency-matrix graph (make-instance 'matrix)))
         (m (matrix-n-rows a))
         (finished)
         (zero 0)
         (one 1))
    (declare (type fixnum one))
    (declare (type fixnum zero))
    (setf nd (make-infinity-matrix nd m m))
    (loop :for i :from 0 :below m :do
       (setf (matrix-ref nd i i) zero))
    (loop :for i :from 0 :below m :do
       (loop :for j :from 0 :below m :do
          (when (= (matrix-ref a i j) one)
            (setf (matrix-ref nd i j) one))))
    (loop :for i :from 2 :to m :unless finished :do
       (setf a-power (matrix-product a a-power))
       (setf finished t)
       (loop :for j :from 0 :below m :do
          (loop :for k :from 0 :below m :do
             (when (and (infinitep (matrix-ref nd j k) nd)
                        (> (matrix-ref a-power j k) zero))
               (setf (matrix-ref nd j k) (coerce i 'fixnum))
               (setf finished nil)))))
    nd))

(defmethod to-distance-matrix ((graph graph) (nd fast-matrix))
  (let* ((a (to-adjacency-matrix graph (make-instance 'fast-matrix)))
         (a-power (to-adjacency-matrix graph (make-instance 'fast-matrix)))
         (m (matrix-n-rows a))
         (finished)
         (zero 0.0s0)
         (one 1.0s0))
    (declare (type single-float one))
    (declare (type single-float zero))
    (setf nd (make-infinity-matrix nd m m))
    (loop :for i :from 0 :below m :do
       (setf (matrix-ref nd i i) zero))
    (loop :for i :from 0 :below m :do
       (loop :for j :from 0 :below m :do
          (when (= (matrix-ref a i j) one)
            (setf (matrix-ref nd i j) one))))
    (loop :for i :from 2 :to m :unless finished :do
       (setf a-power (matrix-product a a-power))
       (setf finished t)
       (loop :for j :from 0 :below m :do
          (loop :for k :from 0 :below m :do
             (when (and (infinitep (matrix-ref nd j k) nd)
                        (> (matrix-ref a-power j k) zero))
               (setf (matrix-ref nd j k) (coerce i 'single-float))
               (setf finished nil)))))
    nd))

(defgeneric distance-from-to (graph nd from to)
  (:documentation "Returns the number of edges in graph GRAPH from
  node FROM to node TO, given the distance matrix ND."))

(defmethod distance-from-to ((graph graph) (nd matrix) from to)
  (let ((node-index-hash (make-hash-table))
        (counter -1))
    (mapc (lambda (node) (setf (gethash node node-index-hash) (incf counter)))
          (nodes graph))
    (matrix-ref nd
                (gethash from node-index-hash)
                (gethash to node-index-hash))))

;; Peirce's relational properties

(defun reflexivep (graph matrix)
  :documentation "Returns t if GRAPH is reflexive, nil otherwise."
  (let ((a (to-adjacency-matrix graph matrix))
        (result))
    (loop :for j :from 0 :below (matrix-n-rows a) :unless result :do
       (setf result (not (= 1 (matrix-ref a j j)))))
    (not result)))

(defun irreflexivep (graph matrix)
  :documentation "Returns t if GRAPH is irreflexive, nil otherwise."
  (let ((a (to-adjacency-matrix graph matrix))
        (result))
    (loop :for j :from 0 :below (matrix-n-rows a) :unless result :do
       (setf result (not (= 0 (matrix-ref a j j)))))
    (not result)))

(defun symmetricp (graph matrix)
  :documentation "Returns t if GRAPH is symmetric, nil otherwise."
  (let* ((a (to-adjacency-matrix graph matrix))
         (at (matrix-transpose a)))
    (if (matrix-entries-different-p a at) nil t)))

(defun asymmetricp (graph matrix)
  :documentation "Returns t if GRAPH is asymmetric, nil otherwise."
  (let* ((a (to-adjacency-matrix graph matrix))
         (at (matrix-transpose a))
         (result))
    (loop :for j :from 0 :below (matrix-n-rows a) :unless result :do
       (loop :for k :from 0 :below (matrix-n-rows a) :unless result :do
          (setf result (and (not (eq j k))
                            (eq (matrix-ref a j k) 1)
                            (eq (matrix-ref at j k) 1)))))
    (not result)))

(defun transitivep (graph matrix)
  :documentation "Returns t if GRAPH is transitive, nil otherwise."
  (let* ((a (to-adjacency-matrix graph matrix))
         (a2 (matrix-product a a))
         (two-path)
         (no-match))
    (loop :for j :from 0 :below (matrix-n-rows a) :do
       (loop :for k :from 0 :below (matrix-n-rows a) :do
          (and (not two-path) (eq (matrix-ref a2 j k) 1)
               (setf two-path j))
          (and (not no-match) (eq (matrix-ref a2 j k) 1)
               (eq (matrix-ref a j k) 0)
               (setf no-match j))))
    (and two-path (not no-match))))

(defun intransitivep (graph matrix)
  :documentation "Returns t if GRAPH is intransitive, nil otherwise."
  (let* ((a (to-adjacency-matrix graph matrix))
         (a2 (matrix-product a a))
         (two-path)
         (match))
    (loop :for j :from 0 :below (matrix-n-rows a) :do
       (loop :for k :from 0 :below (matrix-n-rows a) :do
          (and (not two-path) (eq (matrix-ref a2 j k) 1)
               (setf two-path j))
          (and (not match) (not (eq j k))
               (eq (matrix-ref a2 j k) 1)
               (eq (matrix-ref a j k) 1)
               (setf match j))))
    (and two-path (not match))))

(defun completep (graph matrix)
  :documentation "Returns t if GRAPH is complete, nil otherwise."
  (let* ((a (to-adjacency-matrix graph matrix))
         (at (matrix-transpose a))
         (result))
    (loop :for j :from 0 :below (matrix-n-rows a) :unless result :do
       (loop :for k :from 0 :below (matrix-n-rows a)
          :unless (or result (eq j k)) :do
          (setf result (and (eq (graph/matrix::matrix-ref a j k) 0)
                            (eq (graph/matrix::matrix-ref at j k) 0)))))
    (not result)))

(defun relational-structure (graph matrix)
  :documentation "Returns a string with the name of a relational
  structure whose axiom system GRAPH satisfies, or nil if no
  relational structure axiom system is satisfied."
  (let ((rnobar (reflexivep graph matrix))
        (rbar (irreflexivep graph matrix))
        (snobar (symmetricp graph matrix))
        (sbar (asymmetricp graph matrix))
        (tnobar (transitivep graph matrix))
        (tbar (intransitivep graph matrix))
        (cnobar (completep graph matrix)))
    (or
     (when (and rbar
                (not (or snobar sbar tnobar tbar cnobar))) "digraph")
     (when (and rbar snobar
                (not (or tnobar tbar cnobar))) "graph")
     (when (and rbar sbar
                (not (or tnobar tbar cnobar))) "oriented graph")
     (when (and rnobar snobar
                (not tnobar)) "similarity relation")
     (when (and rnobar snobar tnobar
                (not cnobar)) "equivalence relation")
     (when (and rbar sbar tnobar
                (not cnobar)) "partial order")
     (when (and rbar sbar tnobar cnobar) "complete order")
     (when (and rbar sbar cnobar
                (not (or tnobar tbar))) "tournament")
     (when (and rbar snobar tnobar
                (not cnobar)) "parity relation")
     (when (and rbar sbar tbar
                (not cnobar)) "antiequivalence relation")
     (when (and rnobar sbar tbar) "antiparity relation"))))