817720023

Build # 817720023

Build Type

github

Committed by GitHub

Commit Message

Merge 0e05446eb into 93e37c4d7

Pull Request Pull Request #87: Add Aggregators (WIP)!

Run Details

1725 of 1772 branches covered (97.35%)

Branch coverage included in aggregate %.

597 of 632 new or added lines in 13 files covered. (94.46%)

2 existing lines in 2 files now uncovered.

8155 of 8213 relevant lines covered (99.29%)

10.37 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.44

/grblas/agg.py

"""`grblas.agg` is an experimental module for exploring Aggregators.

Aggregators may be used in reduce methods:
    - Matrix.reduce_rows
    - Matrix.reduce_columns
    - Matrix.reduce_scalar
    - Vector.reduce

Aggregators are implemented as recipes that often use monoids or semirings.

Monoid-only aggregators:
    - sum
    - prod
    - all
    - any
    - min
    - max
    - any_value
    - bitwise_all
    - bitwise_any

Semiring aggregators with O(1) dense vector:
    - count
    - count_nonzero
    - count_zero
    - sum_of_squares
    - sum_of_inverses

Semiring aggregators with UnaryOp applied to the final result:
    - hypot
    - logaddexp
    - logaddexp2

Composite aggregators (require multiple aggregation steps):
    - mean
    - ptp, peak-to-peak
    - varp, population variance
    - vars, sample variance
    - stdp, population standard deviation
    - stds, sample standard deviation
    - geometric_mean
    - harmonic_mean

Custom recipes:
    - first
    - last
    - argmin, doesn't work with Matrix.reduce_scalar
    - argmini
    - argminj, doesn't work with Vector.reduce
    - argmax, doesn't work with Matrix.reduce_scalar
    - argmaxi
    - argmaxj, doesn't work with Vector.reduce
"""
import numpy as np
import grblas as gb
from functools import partial


class Aggregator:
    opclass = "Aggregator"

    def __init__(
        self,
        name,
        *,
        initval=False,
        monoid=None,
        semiring=None,
        switch=False,
        semiring2=None,
        finalize=None,
        composite=None,
        custom=None,
    ):
        self.name = name
        self._initval = initval
        self._initdtype = gb.dtypes.lookup_dtype(type(initval))
        self._monoid = monoid
        self._semiring = semiring
        self._semiring2 = semiring2
        self._switch = switch
        self._finalize = finalize
        self._composite = composite
        self._custom = custom

    def __getitem__(self, dtype):
        return TypedAggregator(self, dtype)


class TypedAggregator:
    opclass = "Aggregator"

    def __init__(self, agg, dtype):
        self.name = agg.name
        self._agg = agg
        self._dtype = dtype
        # TODO: should we specialize the monoid or semiring based on dtype?
        if agg._monoid is not None:
            self.return_type = agg._monoid[dtype].return_type
        elif agg._semiring is not None:
            self.return_type = agg._semiring[dtype].return_type
        elif agg._composite is not None:
            self.return_type = "FP64"  # TODO
        elif agg._custom is not None:
            self.return_type = "INT64"  # TODO
        else:
            raise NotImplementedError()

    def _new(self, updater, expr, *, in_composite=False):
        agg = self._agg
        if agg._monoid is not None:
            updater << getattr(expr.args[0], expr.method_name)(agg._monoid)
            if in_composite:
                parent = updater.parent
                if not parent._is_scalar:
                    return parent
                rv = gb.Vector.new(parent.dtype, size=1, name="TODO")
                rv[0] = parent
                return rv
        elif agg._composite is not None:
            results = []
            mask = updater.kwargs.get("mask")
            for cur_agg in agg._composite:
                arg = expr.construct_output(dtype=updater.parent.dtype)
                results.append(cur_agg[self._dtype]._new(arg(mask=mask), expr, in_composite=True))
            final_expr = agg._finalize(*results)
            if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":
                updater << final_expr
            elif expr.cfunc_name.startswith("GrB_Vector_reduce") or expr.cfunc_name.startswith(
                "GrB_Matrix_reduce"
            ):
                final = final_expr.new(name="TODO")
                updater << final.reduce(gb.monoid.any)
            else:
                raise NotImplementedError()
            if in_composite:
                1 / 0
                return updater.parent
        elif agg._custom is not None:
            return agg._custom(self, updater, expr, in_composite=in_composite)
        elif expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":
            # Matrix -> Vector
            A = expr.args[0]
            orig_updater = updater
            if agg._finalize is not None:
                step1 = expr.construct_output(dtype=updater.parent.dtype)
                updater = step1(mask=updater.kwargs.get("mask"))
            if expr.method_name == "reduce_columns":
                A = A.T
            size = A._nrows if agg._switch else A._ncols
            init = gb.Vector.new(agg._initdtype, size=size, name="TODO")
            init[:] = agg._initval  # O(1) dense vector in SuiteSparse 5
            if agg._switch:
                updater << agg._semiring(init @ A.T)
            else:
                updater << agg._semiring(A @ init)
            if agg._finalize is not None:
                orig_updater << agg._finalize(step1)
            if in_composite:
                return orig_updater.parent
        elif expr.cfunc_name.startswith("GrB_Vector_reduce"):
            # Vector -> Scalar
            v = expr.args[0]
            step1 = gb.Vector.new(updater.parent.dtype, size=1)
            init = gb.Matrix.new(agg._initdtype, nrows=v._size, ncols=1, name="TODO")
            init[:, :] = agg._initval  # O(1) dense column vector in SuiteSparse 5
            if agg._switch:
                step1 << agg._semiring(init.T @ v)
            else:
                step1 << agg._semiring[float](v @ init)
            if agg._finalize is not None:
                step1 << agg._finalize(step1)
            if in_composite:
                return step1
            updater << step1.reduce(gb.monoid.any)
        elif expr.cfunc_name.startswith("GrB_Matrix_reduce"):
            # Matrix -> Scalar
            A = expr.args[0]
            # We need to compute in two steps: Matrix -> Vector -> Scalar.
            # This has not been benchmarked or optimized.
            # We may be able to intelligently choose the faster path.
            if agg._switch:
                size1 = A._nrows
                size2 = A._ncols
            else:
                size1 = A._ncols
                size2 = A._nrows
            init1 = gb.Vector.new(agg._initdtype, size=size1, name="TODO")
            init1[:] = agg._initval  # O(1) dense vector in SuiteSparse 5
            step1 = gb.Vector.new(updater.parent.dtype, size=size2)
            if agg._switch:
                step1 << agg._semiring(init1 @ A.T)
            else:
                step1 << agg._semiring(A @ init1)
            init2 = gb.Matrix.new(agg._initdtype, nrows=size2, ncols=1, name="TODO")
            init2[:, :] = agg._initval  # O(1) dense vector in SuiteSparse 5
            step2 = gb.Vector.new(updater.parent.dtype, size=1)
            step2 << agg._semiring2(step1 @ init2)
            if agg._finalize is not None:
                step2 << agg._finalize(step2)
            if in_composite:
                return step2
            updater << step2.reduce(gb.monoid.any)
        else:
            raise NotImplementedError()


# Monoid-only
sum = Aggregator("sum", monoid=gb.monoid.plus)
prod = Aggregator("prod", monoid=gb.monoid.times)
all = Aggregator("all", monoid=gb.monoid.land)
any = Aggregator("any", monoid=gb.monoid.lor)
min = Aggregator("min", monoid=gb.monoid.min)
max = Aggregator("max", monoid=gb.monoid.max)
any_value = Aggregator("any_value", monoid=gb.monoid.any)
bitwise_all = Aggregator("bitwise_all", monoid=gb.monoid.band)
bitwise_any = Aggregator("bitwise_any", monoid=gb.monoid.bor)
# Other monoids: bxnor bxor eq lxnor lxor

# Semiring-only
count = Aggregator("count", semiring=gb.semiring.plus_pair, semiring2=gb.semiring.plus_first)
count_nonzero = Aggregator(
    "count_nonzero", semiring=gb.semiring.plus_isne, semiring2=gb.semiring.plus_first
)
count_zero = Aggregator(
    "count_zero", semiring=gb.semiring.plus_iseq, semiring2=gb.semiring.plus_first
)
sum_of_squares = Aggregator(
    "sum_of_squares", initval=2, semiring=gb.semiring.plus_pow, semiring2=gb.semiring.plus_first
)
sum_of_inverses = Aggregator(
    "sum_of_inverses", initval=-1, semiring=gb.semiring.plus_pow, semiring2=gb.semiring.plus_first
)

# Semiring and finalize
hypot = Aggregator(
    "hypot",
    initval=2,
    semiring=gb.semiring.plus_pow,
    semiring2=gb.semiring.plus_first,
    finalize=gb.unary.sqrt,
)
logaddexp = Aggregator(
    "logaddexp",
    initval=np.e,
    semiring=gb.semiring.plus_pow,
    switch=True,
    semiring2=gb.semiring.plus_first,
    finalize=gb.unary.log,
)
logaddexp2 = Aggregator(
    "logaddexp2",
    initval=2,
    semiring=gb.semiring.plus_pow,
    switch=True,
    semiring2=gb.semiring.plus_first,
    finalize=gb.unary.log2,
)
# Alternatives
# hypot = Aggregator('hypot', monoid=gb.semiring.numpy.hypot)
# logaddexp = Aggregator('logaddexp', monoid=gb.semiring.numpy.logaddexp)
# logaddexp2 = Aggregator('logaddexp2', monoid=gb.semiring.numpy.logaddexp2)


# Composite
def _mean_finalize(c, x):
    return gb.binary.truediv(x & c)


def _ptp_finalize(max, min):
    return gb.binary.minus(max & min)


def _varp_finalize(c, x, x2):
    # <x2> / n - (<x> / n)**2
    x2 << gb.binary.truediv(x2 & c)
    x << gb.binary.truediv(x & c)
    x << gb.binary.pow(x, 2)
    return gb.binary.minus(x2 & x)


def _vars_finalize(c, x, x2):
    # <x2> / (n-1) - <x>**2 / (n * (n-1))
    x << gb.binary.pow(x, 2)
    x << gb.binary.truediv(x & c)
    c << gb.binary.minus(c, 1)
    x << gb.binary.truediv(x & c)
    x2 << gb.binary.truediv(x2 & c)
    return gb.binary.minus(x2 & x)


def _stdp_finalize(c, x, x2):
    x << _varp_finalize(c, x, x2)
    return gb.unary.sqrt(x)


def _stds_finalize(c, x, x2):
    x << _vars_finalize(c, x, x2)
    return gb.unary.sqrt(x)


def _geometric_mean_finalize(c, x):
    c << gb.unary.minv(c)  # XXX: should be float
    return gb.binary.pow(x & c)


def _harmonic_mean_finalize(c, x):
    return gb.binary.truediv(c & x)


mean = Aggregator(
    "mean",
    composite=[count, sum],
    finalize=_mean_finalize,
)
ptp = Aggregator(
    "ptp",
    composite=[max, min],
    finalize=_ptp_finalize,
)
varp = Aggregator(
    "varp",
    composite=[count, sum, sum_of_squares],
    finalize=_varp_finalize,
)
vars = Aggregator(
    "vars",
    composite=[count, sum, sum_of_squares],
    finalize=_vars_finalize,
)
stdp = Aggregator(
    "stdp",
    composite=[count, sum, sum_of_squares],
    finalize=_stdp_finalize,
)
stds = Aggregator(
    "stds",
    composite=[count, sum, sum_of_squares],
    finalize=_stds_finalize,
)
geometric_mean = Aggregator(
    "geometric_mean",
    composite=[count, prod],
    finalize=_geometric_mean_finalize,
)
harmonic_mean = Aggregator(
    "harmonic_mean",
    composite=[count, sum_of_inverses],
    finalize=_harmonic_mean_finalize,
)


# Special recipes
def _argminmaxij(
    agg,
    updater,
    expr,
    *,
    in_composite,
    monoid,
    col_semiring,
    row_semiring,
    allow_vector=True,
):
    if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":
        A = expr.args[0]
        if expr.method_name == "reduce_rows":
            step1 = A.reduce_rows(monoid).new(name="TODO")

            # TODO: use diag
            i, j = step1.to_values()
            D = gb.Matrix.from_values(i, i, j, nrows=A._nrows, ncols=A._nrows, name="TODO")

            masked = gb.semiring.any_eq(D @ A).new(name="TODO")
            masked(mask=masked.V, replace=True) << masked  # Could use select
            init = gb.Vector.new(bool, size=A._ncols, name="TODO")
            init[:] = False  # O(1) dense vector in SuiteSparse 5
            updater << row_semiring(masked @ init)
            if in_composite:
                1 / 0
                return updater.parent
        else:
            step1 = A.reduce_columns(monoid).new(name="TODO")

            # TODO: use diag
            i, j = step1.to_values()
            D = gb.Matrix.from_values(i, i, j, nrows=A._ncols, ncols=A._ncols, name="TODO")

            masked = gb.semiring.any_eq(A @ D).new(name="TODO")
            masked(mask=masked.V, replace=True) << masked  # Could use select
            init = gb.Vector.new(bool, size=A._nrows, name="TODO")
            init[:] = False  # O(1) dense vector in SuiteSparse 5
            updater << col_semiring(init @ masked)
            if in_composite:
                1 / 0
                return updater.parent
    elif expr.cfunc_name.startswith("GrB_Vector_reduce"):
        if not allow_vector:
            # XXX: it would be best to raise upon creation of the expression
            raise ValueError(
                f"Aggregator {agg.name} may not be used with Vector.reduce; "
                f"use {agg.name[:-1]} instead."
            )
        v = expr.args[0]
        step1 = v.reduce(monoid).new(name="TODO")
        masked = gb.binary.eq(v, step1).new(name="TODO")
        masked(mask=masked.V, replace=True) << masked  # Could use select
        init = gb.Matrix.new(bool, nrows=v._size, ncols=1, name="TODO")
        init[:, :] = False  # O(1) dense column vector in SuiteSparse 5
        step2 = gb.Vector.new(updater.parent.dtype, size=1)
        step2 << col_semiring(masked @ init)
        if in_composite:
            1 / 0
            return step2
        updater << step2.reduce(gb.monoid.any)
    elif expr.cfunc_name.startswith("GrB_Matrix_reduce"):
        A = expr.args[0]
        step1 = A.reduce_scalar(monoid).new(name="TODO")

        masked = gb.binary.eq(A, step1).new(name="TODO")
        masked(mask=masked.V, replace=True) << masked  # Could use select
        init = gb.Vector.new(bool, size=A._nrows, name="TODO")
        init[:] = False  # O(1) dense vector in SuiteSparse 5

        # Always choose the one with smallest i
        step2 = gb.semiring.min_secondi(init @ masked).new(name="TODO")
        step3 = step2.reduce(gb.monoid.min)
        if agg.name in {"argmini", "argmaxi"}:
            # We're done!
            if in_composite:
                1 / 0
                # TODO: put scalar in vector
                raise NotImplementedError()
            updater << step3
            return
        i = step3.value
        step4 = gb.Vector.from_values([i], [False], size=A._nrows, name="TODO")
        step5 = col_semiring(step4 @ masked).new(name="TODO")
        step6 = step5.reduce(gb.monoid.min)
        if in_composite:
            1 / 0
            # TODO: put scalar in vector
            raise NotImplementedError()
        updater << step6
    else:
        raise NotImplementedError()


def _argminmax(agg, updater, expr, *, in_composite, monoid):
    if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":
        if expr.method_name == "reduce_rows":
            return _argminmaxij(
                agg,
                updater,
                expr,
                in_composite=in_composite,
                monoid=monoid,
                row_semiring=gb.semiring.min_firstj,
                col_semiring=gb.semiring.min_secondj,
            )
        return _argminmaxij(
            agg,
            updater,
            expr,
            in_composite=in_composite,
            monoid=monoid,
            row_semiring=gb.semiring.min_firsti,
            col_semiring=gb.semiring.min_secondi,
        )
    elif expr.cfunc_name.startswith("GrB_Vector_reduce"):
        return _argminmaxij(
            agg,
            updater,
            expr,
            in_composite=in_composite,
            monoid=monoid,
            row_semiring=gb.semiring.min_firsti,
            col_semiring=gb.semiring.min_secondi,
        )
    elif expr.cfunc_name.startswith("GrB_Matrix_reduce"):
        # XXX: it would be best to raise upon creation of the expression
        raise ValueError(
            f"Aggregator {agg.name} may not be used with Matrix.reduce_scalar; "
            f"use {agg.name}i or {agg.name}j instead."
        )
    else:
        raise NotImplementedError()


# argmini (argmaxi) is the same as argmin (argmax) for vectors
argmini = Aggregator(
    "argmini",
    custom=partial(
        _argminmaxij,
        monoid=gb.monoid.min,
        row_semiring=gb.semiring.min_firsti,
        col_semiring=gb.semiring.min_secondi,
    ),
)
argmaxi = Aggregator(
    "argmaxi",
    custom=partial(
        _argminmaxij,
        monoid=gb.monoid.max,
        row_semiring=gb.semiring.min_firsti,
        col_semiring=gb.semiring.min_secondi,
    ),
)
# argminj, argmaxj don't work on vectors
argminj = Aggregator(
    "argminj",
    custom=partial(
        _argminmaxij,
        monoid=gb.monoid.min,
        row_semiring=gb.semiring.min_firstj,
        col_semiring=gb.semiring.min_secondj,
        allow_vector=False,
    ),
)
argmaxj = Aggregator(
    "argmaxj",
    custom=partial(
        _argminmaxij,
        monoid=gb.monoid.max,
        row_semiring=gb.semiring.min_firstj,
        col_semiring=gb.semiring.min_secondj,
        allow_vector=False,
    ),
)
# These "do the right thing", but don't work with `reduce_scalar`
argmin = Aggregator("argmin", custom=partial(_argminmax, monoid=gb.monoid.min))
argmax = Aggregator("argmax", custom=partial(_argminmax, monoid=gb.monoid.max))


def _first_last(agg, updater, expr, *, in_composite, semiring):
    if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":
        A = expr.args[0]
        if expr.method_name == "reduce_columns":
            A = A.T
        init = gb.Vector.new(bool, size=A._ncols, name="TODO")
        init[:] = False
        step1 = semiring(A @ init).new(name="TODO")
        Is, Js = step1.to_values()
        # TODO: perform these loops in e.g. Cython
        # Populate numpy array
        vals = np.empty(Is.size, dtype=A.dtype.np_type)
        for index, (i, j) in enumerate(zip(Is, Js)):
            vals[index] = A[i, j].value
        # or Vector
        # v = gb.Vector.new(A.dtype, size=A._nrows)
        # for i, j in zip(Is, Js):
        #     v[i] = A[i, j].value
        result = gb.Vector.from_values(Is, vals, name="TODO")
        updater << result
    elif expr.cfunc_name.startswith("GrB_Vector_reduce"):
        v = expr.args[0]
        init = gb.Matrix.new(bool, nrows=v._size, ncols=1, name="TODO")
        init[:, :] = False
        step1 = semiring(v @ init).new(name="TODO")
        index = step1[0].value
        if index is None:
            index = 0
        updater << v[index]
    else:
        raise NotImplementedError()


first = Aggregator("first", custom=partial(_first_last, semiring=gb.op.min_secondi))
last = Aggregator("last", custom=partial(_first_last, semiring=gb.op.max_secondi))

1	"""`grblas.agg` is an experimental module for exploring Aggregators.
2
3	Aggregators may be used in reduce methods:
4	- Matrix.reduce_rows
5	- Matrix.reduce_columns
6	- Matrix.reduce_scalar
7	- Vector.reduce
8
9	Aggregators are implemented as recipes that often use monoids or semirings.
10
11	Monoid-only aggregators:
12	- sum
13	- prod
14	- all
15	- any
16	- min
17	- max
18	- any_value
19	- bitwise_all
20	- bitwise_any
21
22	Semiring aggregators with O(1) dense vector:
23	- count
24	- count_nonzero
25	- count_zero
26	- sum_of_squares
27	- sum_of_inverses
28
29	Semiring aggregators with UnaryOp applied to the final result:
30	- hypot
31	- logaddexp
32	- logaddexp2
33
34	Composite aggregators (require multiple aggregation steps):
35	- mean
36	- ptp, peak-to-peak
37	- varp, population variance
38	- vars, sample variance
39	- stdp, population standard deviation
40	- stds, sample standard deviation
41	- geometric_mean
42	- harmonic_mean
43
44	Custom recipes:
45	- first
46	- last
47	- argmin, doesn't work with Matrix.reduce_scalar
48	- argmini
49	- argminj, doesn't work with Vector.reduce
50	- argmax, doesn't work with Matrix.reduce_scalar
51	- argmaxi
52	- argmaxj, doesn't work with Vector.reduce
53	"""
54	import numpy as np	11✔
55	import grblas as gb	11✔
56	from functools import partial	11✔
57
58
59	class Aggregator:	11✔
60	opclass = "Aggregator"	11✔
61
62	def __init__(	11✔
63	self,
64	name,
65	*,
66	initval=False,
67	monoid=None,
68	semiring=None,
69	switch=False,
70	semiring2=None,
71	finalize=None,
72	composite=None,
73	custom=None,
74	):
75	self.name = name	11✔
76	self._initval = initval	11✔
77	self._initdtype = gb.dtypes.lookup_dtype(type(initval))	11✔
78	self._monoid = monoid	11✔
79	self._semiring = semiring	11✔
80	self._semiring2 = semiring2	11✔
81	self._switch = switch	11✔
82	self._finalize = finalize	11✔
83	self._composite = composite	11✔
84	self._custom = custom	11✔
85
86	def __getitem__(self, dtype):	11✔
87	return TypedAggregator(self, dtype)	11✔
88
89
90	class TypedAggregator:	11✔
91	opclass = "Aggregator"	11✔
92
93	def __init__(self, agg, dtype):	11✔
94	self.name = agg.name	11✔
95	self._agg = agg	11✔
96	self._dtype = dtype	11✔
97	# TODO: should we specialize the monoid or semiring based on dtype?
98	if agg._monoid is not None:	11✔
99	self.return_type = agg._monoid[dtype].return_type	11✔
100	elif agg._semiring is not None:	11✔
101	self.return_type = agg._semiring[dtype].return_type	11✔
102	elif agg._composite is not None:	11✔
103	self.return_type = "FP64" # TODO	11✔
104	elif agg._custom is not None:	11✔
105	self.return_type = "INT64" # TODO	11✔
106	else:
107	raise NotImplementedError()
108
109	def _new(self, updater, expr, *, in_composite=False):	11✔
110	agg = self._agg	11✔
111	if agg._monoid is not None:	11✔
112	updater << getattr(expr.args[0], expr.method_name)(agg._monoid)	11✔
113	if in_composite:	11✔
114	parent = updater.parent	11✔
115	if not parent._is_scalar:	11✔
116	return parent	11✔
117	rv = gb.Vector.new(parent.dtype, size=1, name="TODO")	11✔
118	rv[0] = parent	11✔
119	return rv	11✔
120	elif agg._composite is not None:	11✔
121	results = []	11✔
122	mask = updater.kwargs.get("mask")	11✔
123	for cur_agg in agg._composite:	11✔
124	arg = expr.construct_output(dtype=updater.parent.dtype)	11✔
125	results.append(cur_agg[self._dtype]._new(arg(mask=mask), expr, in_composite=True))	11✔
126	final_expr = agg._finalize(*results)	11✔
127	if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":	11✔
128	updater << final_expr	11✔
129	elif expr.cfunc_name.startswith("GrB_Vector_reduce") or expr.cfunc_name.startswith(	11✔
130	"GrB_Matrix_reduce"
131	):
132	final = final_expr.new(name="TODO")	11✔
133	updater << final.reduce(gb.monoid.any)	11✔
134	else:
135	raise NotImplementedError()
136	if in_composite:	11!
NEW 137	1 / 0	×
NEW 138	return updater.parent	×
139	elif agg._custom is not None:	11✔
140	return agg._custom(self, updater, expr, in_composite=in_composite)	11✔
141	elif expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":	11✔
142	# Matrix -> Vector
143	A = expr.args[0]	11✔
144	orig_updater = updater	11✔
145	if agg._finalize is not None:	11✔
146	step1 = expr.construct_output(dtype=updater.parent.dtype)	11✔
147	updater = step1(mask=updater.kwargs.get("mask"))	11✔
148	if expr.method_name == "reduce_columns":	11✔
149	A = A.T	11✔
150	size = A._nrows if agg._switch else A._ncols	11✔
151	init = gb.Vector.new(agg._initdtype, size=size, name="TODO")	11✔
152	init[:] = agg._initval # O(1) dense vector in SuiteSparse 5	11✔
153	if agg._switch:	11✔
154	updater << agg._semiring(init @ A.T)	11✔
155	else:
156	updater << agg._semiring(A @ init)	11✔
157	if agg._finalize is not None:	11✔
158	orig_updater << agg._finalize(step1)	11✔
159	if in_composite:	11✔
160	return orig_updater.parent	11✔
161	elif expr.cfunc_name.startswith("GrB_Vector_reduce"):	11✔
162	# Vector -> Scalar
163	v = expr.args[0]	11✔
164	step1 = gb.Vector.new(updater.parent.dtype, size=1)	11✔
165	init = gb.Matrix.new(agg._initdtype, nrows=v._size, ncols=1, name="TODO")	11✔
166	init[:, :] = agg._initval # O(1) dense column vector in SuiteSparse 5	11✔
167	if agg._switch:	11✔
168	step1 << agg._semiring(init.T @ v)	11✔
169	else:
170	step1 << agg._semiring[float](v @ init)	11✔
171	if agg._finalize is not None:	11✔
172	step1 << agg._finalize(step1)	11✔
173	if in_composite:	11✔
174	return step1	11✔
175	updater << step1.reduce(gb.monoid.any)	11✔
176	elif expr.cfunc_name.startswith("GrB_Matrix_reduce"):	11✔
177	# Matrix -> Scalar
178	A = expr.args[0]	11✔
179	# We need to compute in two steps: Matrix -> Vector -> Scalar.
180	# This has not been benchmarked or optimized.
181	# We may be able to intelligently choose the faster path.
182	if agg._switch:	11✔
183	size1 = A._nrows	11✔
184	size2 = A._ncols	11✔
185	else:
186	size1 = A._ncols	11✔
187	size2 = A._nrows	11✔
188	init1 = gb.Vector.new(agg._initdtype, size=size1, name="TODO")	11✔
189	init1[:] = agg._initval # O(1) dense vector in SuiteSparse 5	11✔
190	step1 = gb.Vector.new(updater.parent.dtype, size=size2)	11✔
191	if agg._switch:	11✔
192	step1 << agg._semiring(init1 @ A.T)	11✔
193	else:
194	step1 << agg._semiring(A @ init1)	11✔
195	init2 = gb.Matrix.new(agg._initdtype, nrows=size2, ncols=1, name="TODO")	11✔
196	init2[:, :] = agg._initval # O(1) dense vector in SuiteSparse 5	11✔
197	step2 = gb.Vector.new(updater.parent.dtype, size=1)	11✔
198	step2 << agg._semiring2(step1 @ init2)	11✔
199	if agg._finalize is not None:	11✔
200	step2 << agg._finalize(step2)	11✔
201	if in_composite:	11✔
202	return step2	11✔
203	updater << step2.reduce(gb.monoid.any)	11✔
204	else:
205	raise NotImplementedError()
206
207
208	# Monoid-only
209	sum = Aggregator("sum", monoid=gb.monoid.plus)	11✔
210	prod = Aggregator("prod", monoid=gb.monoid.times)	11✔
211	all = Aggregator("all", monoid=gb.monoid.land)	11✔
212	any = Aggregator("any", monoid=gb.monoid.lor)	11✔
213	min = Aggregator("min", monoid=gb.monoid.min)	11✔
214	max = Aggregator("max", monoid=gb.monoid.max)	11✔
215	any_value = Aggregator("any_value", monoid=gb.monoid.any)	11✔
216	bitwise_all = Aggregator("bitwise_all", monoid=gb.monoid.band)	11✔
217	bitwise_any = Aggregator("bitwise_any", monoid=gb.monoid.bor)	11✔
218	# Other monoids: bxnor bxor eq lxnor lxor
219
220	# Semiring-only
221	count = Aggregator("count", semiring=gb.semiring.plus_pair, semiring2=gb.semiring.plus_first)	11✔
222	count_nonzero = Aggregator(	11✔
223	"count_nonzero", semiring=gb.semiring.plus_isne, semiring2=gb.semiring.plus_first
224	)
225	count_zero = Aggregator(	11✔
226	"count_zero", semiring=gb.semiring.plus_iseq, semiring2=gb.semiring.plus_first
227	)
228	sum_of_squares = Aggregator(	11✔
229	"sum_of_squares", initval=2, semiring=gb.semiring.plus_pow, semiring2=gb.semiring.plus_first
230	)
231	sum_of_inverses = Aggregator(	11✔
232	"sum_of_inverses", initval=-1, semiring=gb.semiring.plus_pow, semiring2=gb.semiring.plus_first
233	)
234
235	# Semiring and finalize
236	hypot = Aggregator(	11✔
237	"hypot",
238	initval=2,
239	semiring=gb.semiring.plus_pow,
240	semiring2=gb.semiring.plus_first,
241	finalize=gb.unary.sqrt,
242	)
243	logaddexp = Aggregator(	11✔
244	"logaddexp",
245	initval=np.e,
246	semiring=gb.semiring.plus_pow,
247	switch=True,
248	semiring2=gb.semiring.plus_first,
249	finalize=gb.unary.log,
250	)
251	logaddexp2 = Aggregator(	11✔
252	"logaddexp2",
253	initval=2,
254	semiring=gb.semiring.plus_pow,
255	switch=True,
256	semiring2=gb.semiring.plus_first,
257	finalize=gb.unary.log2,
258	)
259	# Alternatives
260	# hypot = Aggregator('hypot', monoid=gb.semiring.numpy.hypot)
261	# logaddexp = Aggregator('logaddexp', monoid=gb.semiring.numpy.logaddexp)
262	# logaddexp2 = Aggregator('logaddexp2', monoid=gb.semiring.numpy.logaddexp2)
263
264
265	# Composite
266	def _mean_finalize(c, x):	11✔
267	return gb.binary.truediv(x & c)	11✔
268
269
270	def _ptp_finalize(max, min):	11✔
271	return gb.binary.minus(max & min)	11✔
272
273
274	def _varp_finalize(c, x, x2):	11✔
275	# <x2> / n - (<x> / n)**2
276	x2 << gb.binary.truediv(x2 & c)	11✔
277	x << gb.binary.truediv(x & c)	11✔
278	x << gb.binary.pow(x, 2)	11✔
279	return gb.binary.minus(x2 & x)	11✔
280
281
282	def _vars_finalize(c, x, x2):	11✔
283	# <x2> / (n-1) - <x>*2 / (n (n-1))
284	x << gb.binary.pow(x, 2)	11✔
285	x << gb.binary.truediv(x & c)	11✔
286	c << gb.binary.minus(c, 1)	11✔
287	x << gb.binary.truediv(x & c)	11✔
288	x2 << gb.binary.truediv(x2 & c)	11✔
289	return gb.binary.minus(x2 & x)	11✔
290
291
292	def _stdp_finalize(c, x, x2):	11✔
293	x << _varp_finalize(c, x, x2)	11✔
294	return gb.unary.sqrt(x)	11✔
295
296
297	def _stds_finalize(c, x, x2):	11✔
298	x << _vars_finalize(c, x, x2)	11✔
299	return gb.unary.sqrt(x)	11✔
300
301
302	def _geometric_mean_finalize(c, x):	11✔
303	c << gb.unary.minv(c) # XXX: should be float	11✔
304	return gb.binary.pow(x & c)	11✔
305
306
307	def _harmonic_mean_finalize(c, x):	11✔
308	return gb.binary.truediv(c & x)	11✔
309
310
311	mean = Aggregator(	11✔
312	"mean",
313	composite=[count, sum],
314	finalize=_mean_finalize,
315	)
316	ptp = Aggregator(	11✔
317	"ptp",
318	composite=[max, min],
319	finalize=_ptp_finalize,
320	)
321	varp = Aggregator(	11✔
322	"varp",
323	composite=[count, sum, sum_of_squares],
324	finalize=_varp_finalize,
325	)
326	vars = Aggregator(	11✔
327	"vars",
328	composite=[count, sum, sum_of_squares],
329	finalize=_vars_finalize,
330	)
331	stdp = Aggregator(	11✔
332	"stdp",
333	composite=[count, sum, sum_of_squares],
334	finalize=_stdp_finalize,
335	)
336	stds = Aggregator(	11✔
337	"stds",
338	composite=[count, sum, sum_of_squares],
339	finalize=_stds_finalize,
340	)
341	geometric_mean = Aggregator(	11✔
342	"geometric_mean",
343	composite=[count, prod],
344	finalize=_geometric_mean_finalize,
345	)
346	harmonic_mean = Aggregator(	11✔
347	"harmonic_mean",
348	composite=[count, sum_of_inverses],
349	finalize=_harmonic_mean_finalize,
350	)
351
352
353	# Special recipes
354	def _argminmaxij(	11✔
355	agg,
356	updater,
357	expr,
358	*,
359	in_composite,
360	monoid,
361	col_semiring,
362	row_semiring,
363	allow_vector=True,
364	):
365	if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":	11✔
366	A = expr.args[0]	11✔
367	if expr.method_name == "reduce_rows":	11✔
368	step1 = A.reduce_rows(monoid).new(name="TODO")	11✔
369
370	# TODO: use diag
371	i, j = step1.to_values()	11✔
372	D = gb.Matrix.from_values(i, i, j, nrows=A._nrows, ncols=A._nrows, name="TODO")	11✔
373
374	masked = gb.semiring.any_eq(D @ A).new(name="TODO")	11✔
375	masked(mask=masked.V, replace=True) << masked # Could use select	11✔
376	init = gb.Vector.new(bool, size=A._ncols, name="TODO")	11✔
377	init[:] = False # O(1) dense vector in SuiteSparse 5	11✔
378	updater << row_semiring(masked @ init)	11✔
379	if in_composite:	11!
NEW 380	1 / 0	×
NEW 381	return updater.parent	×
382	else:
383	step1 = A.reduce_columns(monoid).new(name="TODO")	11✔
384
385	# TODO: use diag
386	i, j = step1.to_values()	11✔
387	D = gb.Matrix.from_values(i, i, j, nrows=A._ncols, ncols=A._ncols, name="TODO")	11✔
388
389	masked = gb.semiring.any_eq(A @ D).new(name="TODO")	11✔
390	masked(mask=masked.V, replace=True) << masked # Could use select	11✔
391	init = gb.Vector.new(bool, size=A._nrows, name="TODO")	11✔
392	init[:] = False # O(1) dense vector in SuiteSparse 5	11✔
393	updater << col_semiring(init @ masked)	11✔
394	if in_composite:	11!
NEW 395	1 / 0	×
NEW 396	return updater.parent	×
397	elif expr.cfunc_name.startswith("GrB_Vector_reduce"):	11✔
398	if not allow_vector:	11✔
399	# XXX: it would be best to raise upon creation of the expression
400	raise ValueError(	11✔
401	f"Aggregator {agg.name} may not be used with Vector.reduce; "
402	f"use {agg.name[:-1]} instead."
403	)
404	v = expr.args[0]	11✔
405	step1 = v.reduce(monoid).new(name="TODO")	11✔
406	masked = gb.binary.eq(v, step1).new(name="TODO")	11✔
407	masked(mask=masked.V, replace=True) << masked # Could use select	11✔
408	init = gb.Matrix.new(bool, nrows=v._size, ncols=1, name="TODO")	11✔
409	init[:, :] = False # O(1) dense column vector in SuiteSparse 5	11✔
410	step2 = gb.Vector.new(updater.parent.dtype, size=1)	11✔
411	step2 << col_semiring(masked @ init)	11✔
412	if in_composite:	11!
NEW 413	1 / 0	×
NEW 414	return step2	×
415	updater << step2.reduce(gb.monoid.any)	11✔
416	elif expr.cfunc_name.startswith("GrB_Matrix_reduce"):	11✔
417	A = expr.args[0]	11✔
418	step1 = A.reduce_scalar(monoid).new(name="TODO")	11✔
419
420	masked = gb.binary.eq(A, step1).new(name="TODO")	11✔
421	masked(mask=masked.V, replace=True) << masked # Could use select	11✔
422	init = gb.Vector.new(bool, size=A._nrows, name="TODO")	11✔
423	init[:] = False # O(1) dense vector in SuiteSparse 5	11✔
424
425	# Always choose the one with smallest i
426	step2 = gb.semiring.min_secondi(init @ masked).new(name="TODO")	11✔
427	step3 = step2.reduce(gb.monoid.min)	11✔
428	if agg.name in {"argmini", "argmaxi"}:	11✔
429	# We're done!
430	if in_composite:	11!
NEW 431	1 / 0	×
432	# TODO: put scalar in vector
433	raise NotImplementedError()
434	updater << step3	11✔
435	return	11✔
436	i = step3.value	11✔
437	step4 = gb.Vector.from_values([i], [False], size=A._nrows, name="TODO")	11✔
438	step5 = col_semiring(step4 @ masked).new(name="TODO")	11✔
439	step6 = step5.reduce(gb.monoid.min)	11✔
440	if in_composite:	11!
NEW 441	1 / 0	×
442	# TODO: put scalar in vector
443	raise NotImplementedError()
444	updater << step6	11✔
445	else:
446	raise NotImplementedError()
447
448
449	def _argminmax(agg, updater, expr, *, in_composite, monoid):	11✔
450	if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":	11✔
451	if expr.method_name == "reduce_rows":	11✔
452	return _argminmaxij(	11✔
453	agg,
454	updater,
455	expr,
456	in_composite=in_composite,
457	monoid=monoid,
458	row_semiring=gb.semiring.min_firstj,
459	col_semiring=gb.semiring.min_secondj,
460	)
461	return _argminmaxij(	11✔
462	agg,
463	updater,
464	expr,
465	in_composite=in_composite,
466	monoid=monoid,
467	row_semiring=gb.semiring.min_firsti,
468	col_semiring=gb.semiring.min_secondi,
469	)
470	elif expr.cfunc_name.startswith("GrB_Vector_reduce"):	11✔
471	return _argminmaxij(	11✔
472	agg,
473	updater,
474	expr,
475	in_composite=in_composite,
476	monoid=monoid,
477	row_semiring=gb.semiring.min_firsti,
478	col_semiring=gb.semiring.min_secondi,
479	)
480	elif expr.cfunc_name.startswith("GrB_Matrix_reduce"):	11✔
481	# XXX: it would be best to raise upon creation of the expression
482	raise ValueError(	11✔
483	f"Aggregator {agg.name} may not be used with Matrix.reduce_scalar; "
484	f"use {agg.name}i or {agg.name}j instead."
485	)
486	else:
487	raise NotImplementedError()
488
489
490	# argmini (argmaxi) is the same as argmin (argmax) for vectors
491	argmini = Aggregator(	11✔
492	"argmini",
493	custom=partial(
494	_argminmaxij,
495	monoid=gb.monoid.min,
496	row_semiring=gb.semiring.min_firsti,
497	col_semiring=gb.semiring.min_secondi,
498	),
499	)
500	argmaxi = Aggregator(	11✔
501	"argmaxi",
502	custom=partial(
503	_argminmaxij,
504	monoid=gb.monoid.max,
505	row_semiring=gb.semiring.min_firsti,
506	col_semiring=gb.semiring.min_secondi,
507	),
508	)
509	# argminj, argmaxj don't work on vectors
510	argminj = Aggregator(	11✔
511	"argminj",
512	custom=partial(
513	_argminmaxij,
514	monoid=gb.monoid.min,
515	row_semiring=gb.semiring.min_firstj,
516	col_semiring=gb.semiring.min_secondj,
517	allow_vector=False,
518	),
519	)
520	argmaxj = Aggregator(	11✔
521	"argmaxj",
522	custom=partial(
523	_argminmaxij,
524	monoid=gb.monoid.max,
525	row_semiring=gb.semiring.min_firstj,
526	col_semiring=gb.semiring.min_secondj,
527	allow_vector=False,
528	),
529	)
530	# These "do the right thing", but don't work with `reduce_scalar`
531	argmin = Aggregator("argmin", custom=partial(_argminmax, monoid=gb.monoid.min))	11✔
532	argmax = Aggregator("argmax", custom=partial(_argminmax, monoid=gb.monoid.max))	11✔
533
534
535	def _first_last(agg, updater, expr, *, in_composite, semiring):	11✔
536	if expr.cfunc_name == "GrB_Matrix_reduce_Aggregator":	11✔
537	A = expr.args[0]	11✔
538	if expr.method_name == "reduce_columns":	11✔
539	A = A.T	11✔
540	init = gb.Vector.new(bool, size=A._ncols, name="TODO")	11✔
541	init[:] = False	11✔
542	step1 = semiring(A @ init).new(name="TODO")	11✔
543	Is, Js = step1.to_values()	11✔
544	# TODO: perform these loops in e.g. Cython
545	# Populate numpy array
546	vals = np.empty(Is.size, dtype=A.dtype.np_type)	11✔
547	for index, (i, j) in enumerate(zip(Is, Js)):	11✔
548	vals[index] = A[i, j].value	11✔
549	# or Vector
550	# v = gb.Vector.new(A.dtype, size=A._nrows)
551	# for i, j in zip(Is, Js):
552	# v[i] = A[i, j].value
553	result = gb.Vector.from_values(Is, vals, name="TODO")	11✔
554	updater << result	11✔
555	elif expr.cfunc_name.startswith("GrB_Vector_reduce"):	11✔
556	v = expr.args[0]	11✔
557	init = gb.Matrix.new(bool, nrows=v._size, ncols=1, name="TODO")	11✔
558	init[:, :] = False	11✔
559	step1 = semiring(v @ init).new(name="TODO")	11✔
560	index = step1[0].value	11✔
561	if index is None:	11✔
562	index = 0	11✔
563	updater << v[index]	11✔
564	else:
565	raise NotImplementedError()
566
567
568	first = Aggregator("first", custom=partial(_first_last, semiring=gb.op.min_secondi))	11✔
569	last = Aggregator("last", custom=partial(_first_last, semiring=gb.op.max_secondi))	11✔

metagraph-dev / grblas / 817720023

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous