883408506

Build Type

github

Committed by GitHub

Commit Message

Merge ccd7911d1 into f8902df04

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

Run Details

1847 of 1877 branches covered (98.4%)

Branch coverage included in aggregate %.

421 of 432 new or added lines in 7 files covered. (97.45%)

8 existing lines in 1 file now uncovered.

8449 of 8493 relevant lines covered (99.48%)

20.17 hits per line

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

95.09

/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 __repr__(self):
        return f"{self.name}[{self._dtype}]"

    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)
                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()
                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._ncols
            init = gb.Vector.new(agg._initdtype, size=size)
            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)
            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.
            init1 = gb.Vector.new(agg._initdtype, size=A._ncols)
            init1[:] = agg._initval  # O(1) dense vector in SuiteSparse 5
            step1 = gb.Vector.new(updater.parent.dtype, size=A._nrows)
            if agg._switch:
                step1 << agg._semiring(init1 @ A.T)
            else:
                step1 << agg._semiring(A @ init1)
            init2 = gb.Matrix.new(agg._initdtype, nrows=A._nrows, ncols=1)
            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()

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

            masked = gb.semiring.any_eq(D @ A).new()
            masked(mask=masked.V, replace=True) << masked  # Could use select
            init = gb.Vector.new(bool, size=A._ncols)
            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()

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

            masked = gb.semiring.any_eq(A @ D).new()
            masked(mask=masked.V, replace=True) << masked  # Could use select
            init = gb.Vector.new(bool, size=A._nrows)
            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()
        masked = gb.binary.eq(v, step1).new()
        masked(mask=masked.V, replace=True) << masked  # Could use select
        init = gb.Matrix.new(bool, nrows=v._size, ncols=1)
        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()

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

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

metagraph-dev / grblas / 883408506

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