-
Notifications
You must be signed in to change notification settings - Fork 16
Expand file tree
/
Copy pathmatrix.py
More file actions
3998 lines (3581 loc) · 150 KB
/
Copy pathmatrix.py
File metadata and controls
3998 lines (3581 loc) · 150 KB
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
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import itertools
from collections.abc import Sequence
import numpy as np
from .. import backend, binary, monoid, select, semiring
from ..dtypes import _INDEX, FP64, INT64, lookup_dtype, unify
from ..exceptions import DimensionMismatch, InvalidValue, NoValue, check_status
from . import _supports_udfs, automethods, ffi, lib, utils
from .base import BaseExpression, BaseType, _check_mask, call
from .descriptor import lookup as descriptor_lookup
from .expr import _ALL_INDICES, AmbiguousAssignOrExtract, IndexerResolver, InfixExprBase, Updater
from .mask import Mask, StructuralMask, ValueMask
from .operator import (
UNKNOWN_OPCLASS,
_get_typed_op_from_exprs,
find_opclass,
get_semiring,
get_typed_op,
op_from_string,
)
from .scalar import (
_COMPLETE,
_MATERIALIZE,
Scalar,
ScalarExpression,
ScalarIndexExpr,
_as_scalar,
_scalar_index,
)
from .utils import (
_CArray,
_Pointer,
class_property,
get_order,
ints_to_numpy_buffer,
maybe_integral,
normalize_values,
output_type,
values_to_numpy_buffer,
wrapdoc,
)
from .vector import Vector, VectorExpression, VectorIndexExpr, _isclose_recipe, _select_mask
if backend == "suitesparse":
from .ss.matrix import ss
ffi_new = ffi.new
_CSR_FORMAT = Scalar.from_value(
lib.GrB_CSR_FORMAT, dtype=_INDEX, name="GrB_CSR_FORMAT", is_cscalar=True
)
_CSC_FORMAT = Scalar.from_value(
lib.GrB_CSC_FORMAT, dtype=_INDEX, name="GrB_CSC_FORMAT", is_cscalar=True
)
# COO format is not used yet.
# _COO_FORMAT = Scalar.from_value(
# lib.GrB_COO_FORMAT, dtype=_INDEX, name="GrB_COO_FORMAT", is_cscalar=True
# )
# Custom recipes
def _m_add_v(updater, left, right, op):
full = Vector(right.dtype, left._nrows, name="v_full")
full(**updater.opts)[:] = 0
temp = full.outer(right, binary.second).new(
name="M_temp", mask=updater.kwargs.get("mask"), **updater.opts
)
updater << left.ewise_add(temp, op)
def _m_mult_v(updater, left, right, op):
updater << left.mxm(right.diag(name="M_temp"), get_semiring(monoid.any, op))
def _m_union_m(updater, left, right, left_default, right_default, op):
mask = updater.kwargs.get("mask")
opts = updater.opts
new_left = left.dup(op.type, clear=True)
new_left(mask=mask, **opts) << binary.second(right, left_default)
new_left(mask=mask, **opts) << binary.first(left | new_left)
new_right = right.dup(op.type2, clear=True)
new_right(mask=mask, **opts) << binary.second(left, right_default)
new_right(mask=mask, **opts) << binary.first(right | new_right)
updater << op(new_left & new_right)
def _m_union_v(updater, left, right, left_default, right_default, op):
full = Vector(right.dtype, left._nrows, name="v_full")
full(**updater.opts)[:] = 0
temp = full.outer(right, binary.second).new(
name="M_temp", mask=updater.kwargs.get("mask"), **updater.opts
)
updater << left.ewise_union(temp, op, left_default=left_default, right_default=right_default)
def _reposition(updater, indices, chunk):
updater[indices] = chunk
def _power(updater, A, n, op):
opts = updater.opts
if n == 0:
v = Vector.from_scalar(op.binaryop.monoid.identity, A._nrows, A.dtype, name="v_diag")
updater << v.diag(name="M_diag")
return
if n == 1:
updater << A
return
# Use repeated squaring: compute A^2, A^4, A^8, etc., and combine terms as needed.
# See `numpy.linalg.matrix_power` for a simpler implementation to understand how this works.
# We reuse `result` and `square` outputs, and use `square_expr` so masks can be applied.
result = square = square_expr = None
n, bit = divmod(n, 2)
while True:
if bit != 0:
# Need to multiply `square_expr` or `A` into the result
if square_expr is not None:
# Need to evaluate `square_expr`; either into final result, or into `square`
if n == 0 and result is None:
# Handle `updater << A @ A` without an intermediate value
updater << square_expr
return
if square is None:
# Create `square = A @ A`
square = square_expr.new(name="Squares", **opts)
else:
# Compute `square << square @ square`
square(**opts) << square_expr
square_expr = None
if result is None:
# First time needing the intermediate result!
if square is None:
# Use `A` if possible to avoid unnecessary copying
# We will detect and handle `result is A` below
result = A
else:
# Copy square as intermediate result
result = square.dup(name="Power", **opts)
elif n == 0:
# All done! No more terms to compute
updater << op(result @ square)
return
elif result is A:
# Now we need to create a new matrix for the intermediate result
result = op(result @ square).new(name="Power", **opts)
else:
# Main branch: multiply `square` into `result`
result(**opts) << op(result @ square)
n, bit = divmod(n, 2)
if square_expr is not None:
# We need to perform another squaring, so evaluate current `square_expr` first
if square is None:
# Create `square`
square = square_expr.new(name="Squares", **opts)
else:
# Compute `square`
square << square_expr
if square is None:
# First iteration! Create expression for first square
square_expr = op(A @ A)
else:
# Expression for repeated squaring
square_expr = op(square @ square)
class Matrix(BaseType):
"""Create a new GraphBLAS Sparse Matrix.
Parameters
----------
dtype :
Data type for elements in the Matrix.
nrows : int
Number of rows.
ncols : int
Number of columns.
name : str, optional
Name to give the Matrix. This will be displayed in the ``__repr__``.
"""
__slots__ = "_nrows", "_ncols", "_parent", "ss"
ndim = 2
_is_transposed = False
_name_counter = itertools.count()
__networkx_backend__ = "graphblas"
__networkx_plugin__ = "graphblas"
def __new__(cls, dtype=FP64, nrows=0, ncols=0, *, name=None):
self = object.__new__(cls)
self.dtype = lookup_dtype(dtype)
nrows = _as_scalar(nrows, _INDEX, is_cscalar=True)
ncols = _as_scalar(ncols, _INDEX, is_cscalar=True)
self.name = f"M_{next(Matrix._name_counter)}" if name is None else name
self.gb_obj = ffi_new("GrB_Matrix*")
call("GrB_Matrix_new", [_Pointer(self), self.dtype, nrows, ncols])
self._nrows = nrows.value
self._ncols = ncols.value
self._parent = None
if backend == "suitesparse":
self.ss = ss(self)
return self
@classmethod
def _from_obj(cls, gb_obj, dtype, nrows, ncols, *, parent=None, name=None):
self = object.__new__(cls)
self.gb_obj = gb_obj
self.dtype = dtype
self.name = f"M_{next(Matrix._name_counter)}" if name is None else name
self._nrows = nrows
self._ncols = ncols
self._parent = parent
if backend == "suitesparse":
self.ss = ss(self)
return self
def __del__(self):
parent = getattr(self, "_parent", None)
if parent is not None:
return
gb_obj = getattr(self, "gb_obj", None)
if gb_obj is not None and lib is not None:
# it's difficult/dangerous to record the call, b/c `self.name` may not exist
check_status(lib.GrB_Matrix_free(gb_obj), self)
def _as_vector(self, *, name=None):
"""Cast this Matrix with one column to a Vector.
This is SuiteSparse-specific and may change in the future.
This does not copy the matrix.
"""
if self._ncols != 1:
raise ValueError(
f"Matrix must have a single column (not {self._ncols}) to be cast to a Vector"
)
if backend == "suitesparse":
return Vector._from_obj(
ffi.cast("GrB_Vector*", self.gb_obj),
self.dtype,
self._nrows,
parent=self,
name=f"(GrB_Vector){self.name}" if name is None else name,
)
return self[:, 0].new(name=self.name if name is None else name)
def __repr__(self, mask=None, expr=None):
from .formatting import format_matrix
from .recorder import skip_record
with skip_record:
return format_matrix(self, mask=mask, expr=expr)
def _repr_html_(self, mask=None, collapse=False, expr=None):
if self._parent is not None:
return self._parent._repr_html_(mask=mask, collapse=collapse)
from .formatting import format_matrix_html
from .recorder import skip_record
with skip_record:
return format_matrix_html(self, mask=mask, collapse=collapse, expr=expr)
@property
def _name_html(self):
if self._parent is not None:
return self._parent._name_html
return super()._name_html
def __reduce__(self):
# TODO: we should probably use (or compare to) GraphBLAS serialize methods
if backend == "suitesparse":
pieces = self.ss.export(raw=True)
else:
rows, cols, vals = self.to_coo(sort=False)
pieces = (rows, cols, vals, self.dtype, self._nrows, self._ncols)
return self._deserialize, (pieces, self.name)
@staticmethod
def _deserialize(pieces, name):
if backend == "suitesparse":
return Matrix.ss.import_any(name=name, **pieces)
rows, cols, vals, dtype, nrows, ncols = pieces
return Matrix.from_coo(rows, cols, vals, dtype, nrows=nrows, ncols=ncols, name=name)
@property
def S(self):
"""Create a Mask based on the structure of the Matrix."""
return StructuralMask(self)
@property
def V(self):
"""Create a Mask based on the values in the Matrix (treating each value as truthy)."""
return ValueMask(self)
def __delitem__(self, keys, **opts):
"""Delete a single element, row/column, or submatrix.
Examples
--------
>>> del M[1, 5]
"""
del Updater(self, opts=opts)[keys]
def __getitem__(self, keys):
"""Extract a single element, row/column, or submatrix.
See the `Extract section <../user_guide/operations.html#extract>`__
in the User Guide for more details.
Examples
--------
.. code-block:: python
subM = M[[1, 3, 5], :].new()
"""
resolved_indexes = IndexerResolver(self, keys)
shape = resolved_indexes.shape
if not shape:
return ScalarIndexExpr(self, resolved_indexes)
if len(shape) == 1:
return VectorIndexExpr(self, resolved_indexes, *shape)
nrows, ncols = shape
return MatrixIndexExpr(self, resolved_indexes, nrows, ncols)
def __setitem__(self, keys, expr, **opts):
"""Assign values to a single element, row/column, or submatrix.
See the `Assign section <../user_guide/operations.html#assign>`__
in the User Guide for more details.
Examples
--------
.. code-block:: python
M[0, 0:3] = 17
"""
Updater(self, opts=opts)[keys] = expr
def __contains__(self, index):
"""Indicates whether the (row, col) index has a value present.
Examples
--------
.. code-block:: python
(10, 15) in M
"""
extractor = self[index]
if not extractor._is_scalar:
raise TypeError(
f"Invalid index to Matrix contains: {index!r}. A 2-tuple of ints is expected. "
"Doing `(i, j) in my_matrix` checks whether a value is present at that index."
)
scalar = extractor.new(name="s_contains")
return not scalar._is_empty
def __iter__(self):
"""Iterate over (row, col) indices which are present in the matrix."""
rows, columns, _ = self.to_coo(values=False)
return zip(rows.flat, columns.flat, strict=True)
def __sizeof__(self):
if backend == "suitesparse":
size = ffi_new("size_t*")
check_status(lib.GxB_Matrix_memoryUsage(size, self.gb_obj[0]), self)
return size[0] + object.__sizeof__(self)
raise TypeError("Unable to get size of Matrix with backend: {backend}")
def isequal(self, other, *, check_dtype=False, **opts):
"""Check for exact equality (same size, same structure).
Parameters
----------
other : Matrix
The matrix to compare against
check_dtypes : bool
If True, also checks that dtypes match
Returns
-------
bool
See Also
--------
:meth:`isclose` : For equality check of floating point dtypes
"""
other = self._expect_type(
other, (Matrix, TransposedMatrix), within="isequal", argname="other"
)
if check_dtype and self.dtype != other.dtype:
return False
if self._nrows != other._nrows:
return False
if self._ncols != other._ncols:
return False
if self._nvals != other._nvals:
return False
if check_dtype:
op = binary.eq[self.dtype]
else:
op = get_typed_op(binary.eq, self.dtype, other.dtype, kind="binary")
matches = Matrix(bool, self._nrows, self._ncols, name="M_isequal")
matches(**opts) << self.ewise_mult(other, op)
# ewise_mult performs intersection, so nvals will indicate mismatched empty values
if matches._nvals != self._nvals:
return False
# Check if all results are True
return matches.reduce_scalar(monoid.land, allow_empty=False).new(**opts).value
def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False, **opts):
"""Check for approximate equality (including same size and same structure).
Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``.
Parameters
----------
other : Matrix
The matrix to compare against.
rel_tol : float
Relative tolerance.
abs_tol : float
Absolute tolerance.
check_dtype : bool
If True, also checks that dtypes match
Returns
-------
bool
Whether all values of the Matrix are close to the values in ``other``.
"""
other = self._expect_type(
other, (Matrix, TransposedMatrix), within="isclose", argname="other"
)
if check_dtype and self.dtype != other.dtype:
return False
if self._nrows != other._nrows:
return False
if self._ncols != other._ncols:
return False
if self._nvals != other._nvals:
return False
if not _supports_udfs:
return _isclose_recipe(self, other, rel_tol, abs_tol, **opts)
matches = self.ewise_mult(other, binary.isclose(rel_tol, abs_tol)).new(
bool, name="M_isclose", **opts
)
# ewise_mult performs intersection, so nvals will indicate mismatched empty values
if matches._nvals != self._nvals:
return False
# Check if all results are True
return matches.reduce_scalar(monoid.land, allow_empty=False).new(**opts).value
@property
def nrows(self):
"""Number of rows in the Matrix."""
scalar = _scalar_index("s_nrows")
call("GrB_Matrix_nrows", [_Pointer(scalar), self])
return scalar.gb_obj[0]
@property
def ncols(self):
"""Number of columns in the Matrix."""
scalar = _scalar_index("s_ncols")
call("GrB_Matrix_ncols", [_Pointer(scalar), self])
return scalar.gb_obj[0]
@property
def shape(self):
"""A tuple of ``(nrows, ncols)``."""
return (self._nrows, self._ncols)
@property
def nvals(self):
"""Number of non-empty values in the Matrix."""
scalar = _scalar_index("s_nvals")
call("GrB_Matrix_nvals", [_Pointer(scalar), self])
return scalar.gb_obj[0]
@property
def _nvals(self):
"""Like nvals, but doesn't record calls."""
n = ffi_new("GrB_Index*")
check_status(lib.GrB_Matrix_nvals(n, self.gb_obj[0]), self)
return n[0]
@property
def T(self):
"""Indicates the transpose of the Matrix.
Can be used in the arguments of most operations. It also can be used standalone
as the `Transpose operation <../user_guide/operations.html#transpose>`__.
"""
return TransposedMatrix(self)
def clear(self):
"""In-place operation which clears all values in the Matrix.
After the call, :attr:`nvals` will return 0. The :attr:`shape` will not change.
"""
call("GrB_Matrix_clear", [self])
def resize(self, nrows, ncols):
"""In-place operation which changes the :attr:`shape`.
| Increasing :attr:`nrows` or :attr:`ncols` will expand with empty values.
| Decreasing :attr:`nrows` or :attr:`ncols` will drop existing values above
the new indices.
"""
nrows = _as_scalar(nrows, _INDEX, is_cscalar=True)
ncols = _as_scalar(ncols, _INDEX, is_cscalar=True)
call("GrB_Matrix_resize", [self, nrows, ncols])
self._nrows = nrows.value
self._ncols = ncols.value
def to_coo(self, dtype=None, *, rows=True, columns=True, values=True, sort=True):
"""Extract the indices and values as a 3-tuple of numpy arrays
corresponding to the COO format of the Matrix.
Parameters
----------
dtype :
Requested dtype for the output values array.
rows : bool, default=True
Whether to return rows; will return ``None`` for rows if ``False``
columns :bool, default=True
Whether to return columns; will return ``None`` for columns if ``False``
values : bool, default=True
Whether to return values; will return ``None`` for values if ``False``
sort : bool, default=True
Whether to require sorted indices.
If internally stored rowwise, the sorting will be first by rows, then by column.
If internally stored columnwise, the sorting will be first by column, then by row.
See Also
--------
to_dense
to_edgelist
from_coo
Returns
-------
np.ndarray[dtype=uint64] : Rows
np.ndarray[dtype=uint64] : Columns
np.ndarray : Values
"""
if sort and backend == "suitesparse":
self.wait() # sort in SS
nvals = self._nvals
if rows or backend != "suitesparse":
c_rows = _CArray(size=nvals, name="&rows_array")
else:
c_rows = None
if columns or backend != "suitesparse":
c_columns = _CArray(size=nvals, name="&columns_array")
else:
c_columns = None
if values or backend != "suitesparse":
c_values = _CArray(size=nvals, dtype=self.dtype, name="&values_array")
else:
c_values = None
scalar = _scalar_index("s_nvals")
scalar.value = nvals
dtype_name = "UDT" if self.dtype._is_udt else self.dtype.name
call(
f"GrB_Matrix_extractTuples_{dtype_name}",
[c_rows, c_columns, c_values, _Pointer(scalar), self],
)
if values:
c_values = normalize_values(self, c_values.array, dtype)
if sort and backend != "suitesparse":
col = c_columns.array
row = c_rows.array
ind = np.lexsort((col, row)) # sort by rows, then columns
return (
row[ind] if rows else None,
col[ind] if columns else None,
c_values[ind] if values else None,
)
return (
c_rows.array if rows else None,
c_columns.array if columns else None,
c_values if values else None,
)
def to_edgelist(self, dtype=None, *, values=True, sort=True):
"""Extract the indices and values as a 2-tuple of numpy arrays.
This calls ``to_coo`` then transforms the data into an edgelist.
Parameters
----------
dtype :
Requested dtype for the output values array.
values : bool, default=True
Whether to return values; will return ``None`` for values if ``False``
sort : bool, default=True
Whether to require sorted indices.
If internally stored rowwise, the sorting will be first by rows, then by column.
If internally stored columnwise, the sorting will be first by column, then by row.
See Also
--------
to_coo
to_dense
from_edgelist
Returns
-------
np.ndarray[dtype=uint64] : Edgelist
np.ndarray : Values
"""
rows, columns, values = self.to_coo(dtype, values=values, sort=sort)
return (np.column_stack([rows, columns]), values)
def build(self, rows, columns, values, *, dup_op=None, clear=False, nrows=None, ncols=None):
"""Rarely used method to insert values into an existing Matrix.
The typical use case is to create a new Matrix and insert values
at the same time using :meth:`from_coo`.
All the arguments are used identically in :meth:`from_coo`, except for ``clear``, which
indicates whether to clear the Matrix prior to adding the new values.
"""
# TODO: accept `dtype` keyword to match the dtype of `values`?
rows = ints_to_numpy_buffer(rows, np.uint64, name="row indices")
columns = ints_to_numpy_buffer(columns, np.uint64, name="column indices")
values, _dtype = values_to_numpy_buffer(values, self.dtype)
n = values.shape[0]
if rows.size != n or columns.size != n:
raise ValueError(
"`rows` and `columns` and `values` lengths must match: "
f"{rows.size}, {columns.size}, {values.size}"
)
if clear:
self.clear()
if nrows is not None or ncols is not None:
if nrows is None:
nrows = self._nrows
if ncols is None:
ncols = self._ncols
self.resize(nrows, ncols)
if n == 0:
return
dup_op_given = dup_op is not None
if not dup_op_given:
if not self.dtype._is_udt:
dup_op = binary.plus
elif backend != "suitesparse":
dup_op = binary.any
# SS:SuiteSparse-specific: we use NULL for dup_op
if dup_op is not None:
dup_op = get_typed_op(dup_op, self.dtype, kind="binary")
if dup_op.opclass == "Monoid":
dup_op = dup_op.binaryop
else:
self._expect_op(dup_op, "BinaryOp", within="build", argname="dup_op")
rows = _CArray(rows)
columns = _CArray(columns)
values = _CArray(values, self.dtype)
dtype_name = "UDT" if self.dtype._is_udt else self.dtype.name
call(
f"GrB_Matrix_build_{dtype_name}",
[self, rows, columns, values, _as_scalar(n, _INDEX, is_cscalar=True), dup_op],
)
# Check for duplicates when dup_op was not provided
if not dup_op_given and self._nvals < n:
raise ValueError("Duplicate indices found, must provide `dup_op` BinaryOp")
def dup(self, dtype=None, *, clear=False, mask=None, name=None, **opts):
"""Create a duplicate of the Matrix.
This is a full copy, not a view on the original.
Parameters
----------
dtype :
Data type of the new Matrix. Normal typecasting rules apply.
clear : bool, default=False
If True, the returned Matrix will be empty.
mask : Mask, optional
Mask controlling which elements of the original to
include in the copy.
name : str, optional
Name to give the Matrix.
Returns
-------
Matrix
"""
if dtype is not None or mask is not None or clear:
if dtype is None:
dtype = self.dtype
rv = Matrix(dtype, nrows=self._nrows, ncols=self._ncols, name=name)
if not clear:
rv(mask=mask, **opts)[...] = self
else:
if opts:
# Ignore opts for now
desc = descriptor_lookup(**opts) # noqa: F841 (keep desc in scope for context)
new_mat = ffi_new("GrB_Matrix*")
rv = Matrix._from_obj(new_mat, self.dtype, self._nrows, self._ncols, name=name)
call("GrB_Matrix_dup", [_Pointer(rv), self])
return rv
def diag(self, k=0, dtype=None, *, name=None, **opts):
"""Return a Vector built from the diagonal values of the Matrix.
Parameters
----------
k : int
Off-diagonal offset.
dtype :
Data type of the new Vector. Normal typecasting rules apply.
name : str, optional
Name to give the new Vector.
Returns
-------
:class:`~graphblas.Vector`
"""
if backend == "suitesparse":
from ..ss._core import diag
return diag(self, k=k, dtype=dtype, name=name)
# GraphBLAS spec could use GrB_Vector_diag
if dtype is None:
dtype = self.dtype
if type(k) is not Scalar:
k = Scalar.from_value(k, INT64, is_cscalar=True, name="")
D = select.diag(self, k).new(dtype, name="Diag_temp", **opts)
k = k.value
if k < 0:
size = max(0, min(self._nrows + k, self._ncols))
else:
size = max(0, min(self._ncols - k, self._nrows))
rv = Vector(dtype, size=size, name=name)
if k == 0:
rv(**opts) << D.reduce_rowwise(monoid.any)
else:
d = D.reduce_rowwise(monoid.any).new(name="diag_temp", **opts)
if k < 0:
rv(**opts) << d[d._size - rv._size :]
else:
rv(**opts) << d[: rv._size]
return rv
def wait(self, how="materialize"):
"""Wait for a computation to complete or establish a "happens-before" relation.
Parameters
----------
how : {"materialize", "complete"}
"materialize" fully computes an object.
"complete" establishes a "happens-before" relation useful with multi-threading.
See GraphBLAS documentation for more details.
In `non-blocking mode <../user_guide/init.html#graphblas-modes>`__,
the computations may be delayed and not yet safe to use by multiple threads.
Use wait to force completion of the Matrix.
Has no effect in `blocking mode <../user_guide/init.html#graphblas-modes>`__.
"""
how = how.lower()
if how == "materialize":
mode = _MATERIALIZE
elif how == "complete":
mode = _COMPLETE
else:
raise ValueError(f'`how` argument must be "materialize" or "complete"; got {how!r}')
call("GrB_Matrix_wait", [self, mode])
return self
def get(self, row, col, default=None):
"""Get an element at (``row``, ``col``) indices as a Python scalar.
Parameters
----------
row : int
Row index
col : int
Column index
default :
Value returned if no element exists at (row, col)
Returns
-------
Python scalar
"""
expr = self[row, col]
if expr._is_scalar:
rv = expr.new().value
return default if rv is None else rv
raise ValueError(
"Bad row, col arguments in Matrix.get(...). "
"Indices should get a single element, which will be extracted as a Python scalar."
)
@classmethod
def from_coo(
cls,
rows,
columns,
values=1.0,
dtype=None,
*,
nrows=None,
ncols=None,
dup_op=None,
name=None,
):
"""Create a new Matrix from row and column indices and values.
Parameters
----------
rows : list or np.ndarray
Row indices.
columns : list or np.ndarray
Column indices.
values : list or np.ndarray or scalar, default 1.0
List of values. If a scalar is provided, all values will be set to this single value.
dtype :
Data type of the Matrix. If not provided, the values will be inspected
to choose an appropriate dtype.
nrows : int, optional
Number of rows in the Matrix. If not provided, ``nrows`` is computed
from the maximum row index found in ``rows``.
ncols : int, optional
Number of columns in the Matrix. If not provided, ``ncols`` is computed
from the maximum column index found in ``columns``.
dup_op : :class:`~graphblas.core.operator.BinaryOp`, optional
Function used to combine values if duplicate indices are found.
Leaving ``dup_op=None`` will raise an error if duplicates are found.
name : str, optional
Name to give the Matrix.
See Also
--------
from_dense
from_edgelist
to_coo
Returns
-------
Matrix
"""
rows = ints_to_numpy_buffer(rows, np.uint64, name="row indices")
columns = ints_to_numpy_buffer(columns, np.uint64, name="column indices")
values, dtype = values_to_numpy_buffer(values, dtype, subarray_after=1)
# Compute nrows and ncols if not provided
if nrows is None:
if rows.size == 0:
raise ValueError("No row indices provided. Unable to infer nrows.")
nrows = int(rows.max()) + 1
if ncols is None:
if columns.size == 0:
raise ValueError("No column indices provided. Unable to infer ncols.")
ncols = int(columns.max()) + 1
# Create the new matrix
C = cls(dtype, nrows, ncols, name=name)
if values.ndim == 0:
if dup_op is not None:
raise ValueError(
"dup_op must be None if values is a scalar so that all "
"values can be identical. Duplicate indices will be ignored."
)
if backend == "suitesparse":
C.ss.build_scalar(rows, columns, values.tolist())
else:
C.build(rows, columns, np.broadcast_to(values, rows.size), dup_op=binary.any)
else:
# Add the data
# This needs to be the original data to get proper error messages
C.build(rows, columns, values, dup_op=dup_op)
return C
@classmethod
def from_edgelist(
cls,
edgelist,
values=None,
dtype=None,
*,
nrows=None,
ncols=None,
dup_op=None,
name=None,
):
"""Create a new Matrix from edgelist of (row, col) pairs or (row, col, value) triples.
This transforms the data and calls ``Matrix.from_coo``.
Parameters
----------
edgelist : list or np.ndarray or iterable
A sequence of ``(row, column)`` pairs or ``(row, column, value)`` triples.
NumPy edgelist only supports ``(row, column)``; values must be passed separately.
values : list or np.ndarray or scalar, optional
List of values. If a scalar is provided, all values will be set to this single value.
The default is 1.0 if ``edgelist`` is a sequence of ``(row, column)`` pairs.
dtype :
Data type of the Matrix. If not provided, the values will be inspected
to choose an appropriate dtype.
nrows : int, optional
Number of rows in the Matrix. If not provided, ``nrows`` is computed
from the maximum row index found in the edgelist.
ncols : int, optional
Number of columns in the Matrix. If not provided, ``ncols`` is computed
from the maximum column index found in the edgelist.
dup_op : :class:`~graphblas.core.operator.BinaryOp`, optional
Function used to combine values if duplicate indices are found.
Leaving ``dup_op=None`` will raise an error if duplicates are found.
name : str, optional
Name to give the Matrix.
See Also
--------
from_coo
from_dense
to_edgelist
Returns
-------
Matrix
"""
edgelist_values = None
if isinstance(edgelist, np.ndarray):
if edgelist.ndim != 2:
raise ValueError(
f"edgelist array must have 2 dimensions (nvals x 2); got {edgelist.ndim}"
)
if edgelist.shape[1] == 3:
raise ValueError(
"edgelist as NumPy array only supports ``(row, column)``; "
"values must be passed separately."
)
if edgelist.shape[1] != 2:
raise ValueError(
"Last dimension of edgelist array must be length 2 "
f"(for row and column); got {edgelist.shape[1]}"
)
rows = edgelist[:, 0]
cols = edgelist[:, 1]
else:
unzipped = list(zip(*edgelist, strict=True))
if len(unzipped) == 2:
rows, cols = unzipped
elif len(unzipped) == 3:
rows, cols, edgelist_values = unzipped
elif not unzipped:
# Empty edgelist (nrows and ncols should be given)
rows = cols = unzipped
else:
raise ValueError(
"Each item in the edgelist must have two or three elements "
f"(for row and column index, and maybe values); got {len(unzipped)}"
)
if values is None:
if edgelist_values is None:
values = 1.0
else:
values = edgelist_values
elif edgelist_values is not None:
raise TypeError(
"Too many sources of values: from `edgelist` triples and from `values=` argument"
)
return cls.from_coo(
rows, cols, values, dtype, nrows=nrows, ncols=ncols, dup_op=dup_op, name=name
)
@classmethod
def _from_csx(cls, fmt, indptr, indices, values, dtype, num, check_num, name):
if fmt is _CSR_FORMAT:
indices_name = "column indices"
else:
indices_name = "row indices"
indptr = ints_to_numpy_buffer(indptr, np.uint64, name="index pointers")
indices = ints_to_numpy_buffer(indices, np.uint64, name=indices_name)
values, dtype = values_to_numpy_buffer(values, dtype, subarray_after=1)
if num is None: