+
+
+.. raw:: latex
+
+ \nbsphinxstartnotebook{\scriptsize\noindent\strut
+ \textcolor{gray}{The following section was generated from
+ \sphinxcode{\sphinxupquote{\strut {{ nbname | escape_latex }}}} \dotfill}}
+"""
diff --git a/docs/contributor_guide/index.rst b/docs/contributor_guide/index.rst
new file mode 100644
index 000000000..717a78468
--- /dev/null
+++ b/docs/contributor_guide/index.rst
@@ -0,0 +1,5 @@
+
+Contributor Guide
+=================
+
+Needs work
\ No newline at end of file
diff --git a/docs/getting_started/faq.rst b/docs/getting_started/faq.rst
new file mode 100644
index 000000000..df1b679ab
--- /dev/null
+++ b/docs/getting_started/faq.rst
@@ -0,0 +1,58 @@
+
+.. _faq:
+
+Frequently Asked Questions
+==========================
+
+FAQs
+----
+
+What is the relationship between python-graphblas and pygraphblas?
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+They are different libraries with similar goals. Both wrap SuiteSparse:GraphBLAS and use the
+same underlying ``python-suitesparse-graphblas`` CFFI wrapper, making their underlying objects
+compatible with each other.
+
+The biggest difference is in how the C calls are mapped to equivalent Python functions.
+
+For simple expressions, the two libraries are very similar.
+
+.. code-block:: python
+
+ # pygraphblas
+ A += 1
+
+.. code-block:: python
+
+ # python-graphblas
+ A += 1
+
+For more complex expressions, however, the two libraries diverge significantly.
+`pygraphblas` tends towards the `numpy` style of immediate execution and using
+keyword arguments to affect the output.
+
+.. code-block:: python
+
+ # pygraphblas
+ A.mxm(B.transpose(), mask=A, out=C, accum=FP64.PLUS, semiring=FP64.MIN_PLUS)
+
+`python-graphblas` uses delayed expressions and keeps the output-affecting arguments
+together with the output.
+
+.. code-block:: python
+
+ # python-graphblas
+ C(A, accum=binary.plus) << semiring.min_plus(A @ B.T)
+
+`python-graphblas` also contains additional features, such as the `Recorder` and advanced aggregators.
+
+What is the performance penalty of writing algorithms with python-graphblas vs writing them directly in C?
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+For large graphs, the performance penalty is very small as the bulk of the work
+of operating on Matrices and Vectors is done purely in C. The overhead is only in the calls made
+to C.
+
+For small graphs, the call overhead may become more significant, but smaller graphs usually don't
+take much time to compute, so the extra overhead may not be noticeable at a human scale.
\ No newline at end of file
diff --git a/docs/getting_started/index.rst b/docs/getting_started/index.rst
new file mode 100644
index 000000000..661550803
--- /dev/null
+++ b/docs/getting_started/index.rst
@@ -0,0 +1,53 @@
+
+Getting Started
+===============
+
+Installation
+------------
+
+Using conda:
+
+::
+
+ conda install -c conda-forge python-graphblas
+
+Using pip:
+
+::
+
+ pip install python-graphblas
+
+Whether installing with conda or pip, the underlying package that is imported in Python
+is named ``graphblas``. The convention is to import as:
+
+.. code-block:: python
+
+ >>> import graphblas as gb
+
+Optional Dependencies
++++++++++++++++++++++
+
+The following are not required by ``python-graphblas``, but may be needed for certain functionality
+to work.
+
+ - `pandas `__ -- required for nicer ``__repr__``
+ - `matplotlib `__ -- required for basic plotting of graphs
+ - `scipy `__ -- used in ``io`` module to read/write ``scipy.sparse`` format
+ - `networkx `__ -- used in ``io`` module to interface with networkx graphs
+
+GraphBLAS Fundamentals
+----------------------
+
+For a short introduction to the concepts of graph analytics using linear algebra,
+read the :ref:`primer`.
+
+For more details, the best resource for learning about GraphBLAS is `graphblas.org `_.
+It contains information about the history, core ideas, and the full specification. It also contains links to
+many videos and papers about GraphBLAS, as well as a list of implementations and language wrappers.
+
+.. toctree::
+ :maxdepth: 1
+ :hidden:
+
+ primer
+ faq
diff --git a/docs/getting_started/primer.rst b/docs/getting_started/primer.rst
new file mode 100644
index 000000000..0b19e30ad
--- /dev/null
+++ b/docs/getting_started/primer.rst
@@ -0,0 +1,282 @@
+
+.. _primer:
+
+GraphBLAS Primer
+================
+
+This primer covers:
+
+ - :ref:`what-is-a-graph`
+ - :ref:`common-graph-storage`
+ - :ref:`pagerank-using-linear-algebra`
+ - :ref:`introduction-to-semirings`
+ - :ref:`sssp-in-python-graphblas`
+ - :ref:`summary`
+
+.. _what-is-a-graph:
+
+What is a graph?
+----------------
+
+A graph is a set of nodes (also called vertices) that are connected with edges.
+
+When the edges have an associated direction (i.e. arrow), the graph is *directed*.
+
+The following graph is a task graph, meaning subsequent tasks can only be started
+when all previous dependencies are completed.
+
+.. image:: ../_static/img/task-graph.svg
+ :alt: Task Graph
+
+When edges do not have a defined direction, the graph is *undirected*.
+
+In the following social network, it is assumed friendships are mutual,
+so direction is unnecessary when indicating friendships as edges.
+
+Notice that the graph is not fully connected. There is no requirement for all nodes in a
+graph to be reachable from any other node.
+
+.. image:: ../_static/img/social-network.svg
+ :alt: Social Network
+
+Edges can have weights associated them. The meaning of the weight depends on the graph.
+For a road network, it might represent the travel time based on the length of road and
+the posted speed limit.
+
+This graph shows the complexity that can occur in a graph. Self-edges have the same node
+as their start and end point.
+
+.. image:: ../_static/img/directed-graph.svg
+ :alt: Directed Graph
+
+
+.. _common-graph-storage:
+
+Common ways to store a graph
+----------------------------
+
+The most common way to store a graph is as an edge list, containing the `(start_node, end_node, weight)`.
+
+The following graph and its edge list are shown below.
+
+.. image:: ../_static/img/undirected-graph.svg
+
+::
+
+ (0, 1, 5.6)
+ (0, 2, 2.3)
+ (0, 3, 4.6)
+ (1, 2, 1.9)
+ (1, 3, 6.2)
+ (2, 3, 3.0)
+ (3, 4, 1.4)
+ (4, 5, 4.4)
+ (4, 6, 1.0)
+ (5, 6, 2.8)
+
+Notice that each edge is only listed once. An edge list would need to indicate whether
+these are directed edges or undirected.
+
+While edge lists are easy to create, they are not ideal for performance. Finding all the
+edges connected to node 3 requires a full scan of the list. For this reason, other formats
+exist which make common operations more performant.
+
+For example, networkx uses a dict of dicts to store graphs, allowing easy query of each node's
+outgoing edges. This requires doubling the memory usage in the case of undirected graphs, but
+makes for faster graph algorithms.
+
+.. code-block:: python
+
+ # networkx-style storage of an undirected graph
+ G = {
+ 0: {1: {'weight': 5.6},
+ 2: {'weight': 2.3},
+ 3: {'weight': 4.6}},
+ 1: {0: {'weight': 5.6},
+ 2: {'weight': 1.9},
+ 3: {'weight': 6.2}},
+ 2: {0: {'weight': 2.3},
+ 1: {'weight': 1.9},
+ 3: {'weight': 3.0}},
+ 3: {0: {'weight': 4.6},
+ 1: {'weight': 6.2},
+ 2: {'weight': 3.0},
+ 4: {'weight': 1.4}},
+ 4: {3: {'weight': 1.4},
+ 5: {'weight': 4.4},
+ 6: {'weight': 1.0}},
+ 5: {4: {'weight': 4.4},
+ 6: {'weight': 2.8}},
+ 6: {4: {'weight': 1.0},
+ 5: {'weight': 2.8}}
+ }
+
+An alternative way to store a graph is as an adjacency matrix. Each node becomes both a row
+and a column in the matrix. Cells in the matrix are filled in with the edge weight (if it exists).
+This naturally creates a square sparse matrix.
+
+.. image:: ../_static/img/adj-graph.png
+ :alt: Sparse Adjacency Graph
+
+Reading along a row, the outgoing edges are shown, while reading down a column show the
+incoming edges to a node.
+
+In the case of an undirected graph, the graph is symmetric, so the property may not be apparent
+for this example.
+
+While the adjacency matrix is a nice format to work with conceptually, it would be very inefficient
+to store the full matrix as a 2-dimensional dense array. Instead, formats exist to efficiently store
+only the present values.
+
+Various formats exist for efficiently storing only the present values. Examples include
+Compressed Sparse Row (CSR) as well as hypersparse CSR format (also called double-compressed
+sparse row). These are very good at iterating over the sparse matrix in a row-wise manner.
+Similar formats exist which are column-oriented.
+
+.. _pagerank-using-linear-algebra:
+
+Pagerank using Linear Algebra
+-----------------------------
+
+Once we enter the realm of sparse matrices for storing graphs, the natural question is whether
+the field of linear algebra can be used to do anything useful with the graph representation.
+
+Let's look at a motivating example with Pagerank, which is one of the most famous graph analysis
+metrics in the area of web search.
+
+Pagerank takes the weight of a given node and distributes it to each of its connecting nodes in
+an iterative manner until the process converges. The final value of each node is the pagerank
+and will generally indicate the importance of that node in the overall graph.
+
+This is a simplified version of the networkx pagerank algorithm. The outermost for-loop is for
+the overall convergence of the algorithm. The two inner for-loops iterate over each edge of each
+node. This is rather slow because of how Python handles for-loops.
+
+.. code-block:: python
+
+ for _ in range(max_iter):
+ xlast = x
+ x = dict.fromkeys(xlast.keys(), 0)
+ for n in x:
+ for _, nbr, wt in W.edges(n, data=weight):
+ x[nbr] += alpha * xlast[n] * wt
+ # check convergence, l1 norm
+ err = sum(abs(x[n] - xlast[n]) for n in x)
+ if err < N * tol:
+ return x
+ raise nx.PowerIterationFailedConvergence(max_iter)
+
+Pagerank has a very natural representation using linear algebra. The inner for-loops are simply
+performing a matrix multiplication with a vector.
+
+Here is the same algorithm using scipy.sparse matrix and vector. It is both simpler to read and
+also much faster, as the matrix multiplication can be dispatched to scipy.
+
+.. code-block:: python
+
+ for _ in range(max_iter):
+ xlast = x
+ x = alpha * (x @ A)
+ # check convergence, l1 norm
+ err = np.absolute(x - xlast).sum()
+ if err < N * tol:
+ return dict(zip(nodelist, map(float, x)))
+ raise nx.PowerIterationFailedConvergence(max_iter)
+
+
+
+.. _introduction-to-semirings:
+
+Introduction to semirings
+-------------------------
+
+With a sparse vector representing the current value of each node and a sparse adjacency matrix
+representing the edges of the graph, a Vector-Matrix multiplication is equivalent to performing
+one step of a Breadth-First Search (BFS).
+
+The value of each node is multiplied by the weight of each edge, and those new values are summed
+together by the receiving node. In this way, the initial node's value spreads out throughout the graph
+first to the one-hop neighbors and then to the two-hop neighbors, etc.
+
+This idea of propagating values to neighbors is fundamental to many graph algorithms. However, rather
+than using the standard matrix multiplication operators of multiplying elements together and adding
+the results, other operators are used in the same manner. These are called **semirings**.
+
+The "standard" way of performing matrix multiplication is using the *plus_times* semiring.
+
+.. image:: ../_static/img/plus-times-semiring.png
+ :alt: sum(k) (a_ik * b_kj)
+
+An alternative semiring is the *min_plus* semiring.
+
+.. image:: ../_static/img/min-plus-semiring.png
+ :alt: min(k) (a_ik + b_kj)
+
+Conceptually, the *min_plus* semiring takes a node's value, *adds* it to the edge weight, and when
+collecting values at the receiving node, it takes the smallest value.
+
+Why is this useful?
+
+Because doing a repeated BFS with the *min_plus* semiring will compute the Single-Source
+Shortest Path (SSSP) algorithm using linear algebra.
+
+.. _sssp-in-python-graphblas:
+
+SSSP in python-graphblas
+------------------------
+
+.. image:: ../_static/img/super-simple.svg
+
+Below is the python-graphblas algorithm to compute Single-Source Shortest Path for this graph, starting from
+node 0.
+
+.. code-block:: python
+
+ from graphblas import op, semiring, Matrix, Vector
+
+ # Create the graph and starting vector
+ start_node = 0
+ G = Matrix.from_values(
+ [0, 0, 1, 1, 2],
+ [1, 2, 2, 3, 3],
+ [2.0, 5.0, 1.5, 4.25, 0.5],
+ nrows=4, ncols=4)
+ v = Vector.from_values([start_node], [0.0], size=4)
+
+ # Compute SSSP
+ while True:
+ w = v.dup()
+
+ # Perform a BFS step using min_plus semiring
+ # Accumulate into v using the `min` function
+ v(op.min) << semiring.min_plus(v @ G)
+
+ # The algorithm is converged once v stops changing
+ if v.isequal(w):
+ break
+
+The result finds the shortest path from node 0 -> node 3 has a path length of 4.0, which
+goes through nodes 1 and 2. It is not the most direct, but it is the shortest.
+
+.. image:: ../_static/img/sssp-result.png
+ :alt: Vector [0.0, 2.0, 3.5, 4.0 ]
+
+Other algorithms exist to identify the parent of each node when walking back to the source node,
+but this algorithm only computes the shortest path length.
+
+.. _summary:
+
+Summary
+-------
+
+This primer has been a short introduction to graphs, representing graphs as sparse adjacency matrices,
+and showing that linear algebra can be used to compute graph algorithms with the expanded concept
+of semirings.
+
+This is a somewhat new field of research, so many academic papers and talks are being given every year.
+`Graphblas.org `_ remains the best source for keeping up-to-date with the latest
+developments in this area.
+
+Many people will benefit from faster graph algorithms written in GraphBLAS, but for those that want
+to play around with the full power of GraphBLAS, there is no substitute for getting your hands dirty
+in code. Read the :ref:`user_guide` to learn more.
diff --git a/docs/index.rst b/docs/index.rst
index dc05d52e8..ba752d558 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,227 +1,78 @@
-.. graphblas documentation master file, created by
- sphinx-quickstart on Tue Mar 15 13:13:40 2022.
- You can adapt this file completely to your liking, but it should at least
- contain the root `toctree` directive.
-python-graphblas - A Python Wrapper for GraphBLAS
-=======================================
+.. module:: python-graphblas
-To install, ``conda install -c conda-forge python-graphblas`` or ``pip install python-graphblas``. This will also install the SuiteSparse ``graphblas`` compiled C library.
+python-graphblas documentation
+==============================
-Currently works with `SuiteSparse:GraphBLAS `_, but the goal is to make it work with all implementations of the GraphBLAS spec.
+:mod:`python-graphblas` is a Pythonic interface to the highly performant
+`SuiteSparse:GraphBLAS `_ library
+for performing graph analytics in the language of linear algebra.
-The approach taken with this library is to follow the C-API specification as closely as possible while making improvements
-allowed with the Python syntax. Because the spec always passes in the output object to be written to, we follow the same,
-which is very different from the way Python normally operates. In fact, many who are familiar with other Python data
-libraries (numpy, pandas, etc) will find it strange to not create new objects for every call.
+Licensed under the `Apache 2.0 license `_.
-At the highest level, the goal is to separate output, mask, and accumulator on the left side of the assignment
-operator ``=`` and put the computation on the right side. Unfortunately, that approach doesn't always work very well
-with how Python handles assignment, so instead we (ab)use the left-shift ``<<`` notation to give the same flavor of
-assignment. This opens up all kinds of nice possibilities.
+.. toctree::
+ :maxdepth: 2
+ :hidden:
-This is an example of how the mapping works:
+ getting_started/index
+ user_guide/index
+ api_reference/index
+ contributor_guide/index
-.. code-block:: c
+.. panels::
+ :card: + intro-card text-center
+ :column: col-lg-6 col-md-6 col-sm-6 col-xs-12 d-flex
- // C call
- GrB_Matrix_mxm(M, mask, GrB_PLUS_INT64, GrB_MIN_PLUS_INT64, A, B, NULL)
+ ---
+ :fa:`play-circle,fa-3x`
-.. code-block:: python
+ Getting Started
+ ^^^^^^^^^^^^^^^
- # Python call
- M(mask.V, accum=binary.plus) << A.mxm(B, semiring.min_plus)
+ .. link-button:: getting_started
+ :text:
+ :classes: stretched-link
-The expression on the right ``A.mxm(B)`` creates a delayed object which does no computation. Once it is used in the
-``<<`` expression with ``M``, the whole thing is translated into the equivalent GraphBLAS call.
+ Exactly what you need to get up and running with
+ python-graphblas. Contains a short introduction
+ to the main ideas behind GraphBLAS.
-Delayed objects also have a ``.new()`` method which can be used to force computation and return a new
-object. This is convenient and often appropriate, but will create many unnecessary objects if used in a loop. It
-also loses the ability to perform accumulation with existing results. For best performance, following the standard
-GraphBLAS approach of (1) creating the object outside the loop and (2) using the object repeatedly within each loop
-is a much better approach, even if it doesn't feel very Pythonic.
+ ---
+ :fa:`code,fa-3x`
-Descriptor flags are set on the appropriate elements to keep logic close to what it affects. Here is the same call
-with descriptor bits set. ``ttcsr`` indicates transpose the first and second matrices, complement the structure of the mask,
-and do a replacement on the output.
+ User Guide
+ ^^^^^^^^^^
-.. code-block:: c
+ .. link-button:: user_guide
+ :text:
+ :classes: stretched-link
- // C call
- GrB_Matrix_mxm(M, mask, GrB_PLUS_INT64, GrB_MIN_PLUS_INT64, A, B, desc.ttcsr)
+ The user guide covers the fundamental ideas and usage
+ patterns of python-graphblas. It also contains an example
+ notebook showing Single-Source Shortest Path.
-.. code-block:: python
+ ---
+ :fa:`cogs,fa-3x`
- # Python call
- M(~mask.S, accum=binary.plus, replace=True) << A.T.mxm(B.T, semiring.min_plus)
+ API Reference
+ ^^^^^^^^^^^^^
-The objects receiving the flag operations (``A.T``, ``~mask``, etc) are also delayed objects. They hold on to the state but
-do no computation, allowing the correct descriptor bits to be set in a single GraphBLAS call.
+ .. link-button:: api_reference
+ :text:
+ :classes: stretched-link
-If no mask or accumulator is used, the call looks like this:
+ The API reference contains details about each class and method.
+ It assumes familiarity with the core ideas and usage patterns.
-.. code-block:: python
+ ---
+ :fa:`plus-circle,fa-3x`
- M << A.mxm(B, semiring.min_plus)
+ Contributor Guide
+ ^^^^^^^^^^^^^^^^^
-The use of ``<<`` to indicate updating is actually just syntactic sugar for a real ``.update()`` method. The above
-expression could be written as:
+ .. link-button:: contributor_guide
+ :text:
+ :classes: stretched-link
-.. code-block:: python
-
- M.update(A.mxm(B, semiring.min_plus))
-
-Operations
-----------
-
-.. code-block:: python
-
- M(mask, accum) << A.mxm(B, semiring) # mxm
- w(mask, accum) << A.mxv(v, semiring) # mxv
- w(mask, accum) << v.vxm(B, semiring) # vxm
- M(mask, accum) << A.ewise_add(B, binaryop) # eWiseAdd
- M(mask, accum) << A.ewise_mult(B, binaryop) # eWiseMult
- M(mask, accum) << A.kronecker(B, binaryop) # kronecker
- M(mask, accum) << A.T # transpose
-
-Extract
--------
-
-.. code-block:: python
-
- M(mask, accum) << A[rows, cols] # rows and cols are a list or a slice
- w(mask, accum) << A[rows, col_index] # extract column
- w(mask, accum) << A[row_index, cols] # extract row
- s = A[row_index, col_index].value # extract single element
-
-Assign
-------
-
-.. code-block:: python
-
- M(mask, accum)[rows, cols] << A # rows and cols are a list or a slice
- M(mask, accum)[rows, col_index] << v # assign column
- M(mask, accum)[row_index, cols] << v # assign row
- M(mask, accum)[rows, cols] << s # assign scalar to many elements
- M[row_index, col_index] << s # assign scalar to single element
- # (mask and accum not allowed)
- del M[row_index, col_index] # remove single element
-
-Apply
------
-
-.. code-block:: python
-
- M(mask, accum) << A.apply(unaryop)
- M(mask, accum) << A.apply(binaryop, left=s) # bind-first
- M(mask, accum) << A.apply(binaryop, right=s) # bind-second
-
-Reduce
-------
-
-.. code-block:: python
-
- v(mask, accum) << A.reduce_rowwise(op) # reduce row-wise
- v(mask, accum) << A.reduce_columnwise(op) # reduce column-wise
- s(accum) << A.reduce_scalar(op)
- s(accum) << v.reduce(op)
-
-Creating new Vectors / Matrices
--------------------------------
-
-.. code-block:: python
-
- A = Matrix.new(dtype, num_rows, num_cols) # new_type
- B = A.dup() # dup
- A = Matrix.from_values([row_indices], [col_indices], [values]) # build
-
-New from delayed
-----------------
-
-Delayed objects can be used to create a new object using ``.new()`` method.
-
-.. code-block:: python
-
- C = A.mxm(B, semiring).new()
-
-Properties
-----------
-
-.. code-block:: python
-
- size = v.size # size
- nrows = M.nrows # nrows
- ncols = M.ncols # ncols
- nvals = M.nvals # nvals
- rindices, cindices, vals = M.to_values() # extractTuples
-
-Initialization
---------------
-
-There is a mechanism to initialize ``graphblas`` with a context prior to use. This allows for setting the backend to
-use as well as the blocking/non-blocking mode. If the context is not initialized, a default initialization will
-be performed automatically.
-
-.. code-block:: python
-
- import graphblas as gb
- # Context initialization must happen before any other imports
- gb.init('suitesparse', blocking=True)
-
- # Now we can import other items from graphblas
- from graphblas import binary, semiring
- from graphblas import Matrix, Vector, Scalar
-
-Performant User Defined Functions
----------------------------------
-
-Python-graphblas requires ``numba`` which enables compiling user-defined Python functions to native machine code for use in GraphBLAS.
-
-Example customized UnaryOp:
-
-.. code-block:: python
-
- from graphblas import unary
- from graphblas.operator import UnaryOp
-
- def force_odd_func(x):
- if x % 2 == 0:
- return x + 1
- return x
-
- UnaryOp.register_new('force_odd', force_odd_func)
-
- v = Vector.from_values([0, 1, 3], [1, 2, 3])
- w = v.apply(unary.force_odd).new()
- w # indexes=[0, 1, 3], values=[1, 3, 3]
-
-Similar methods exist for BinaryOp, Monoid, and Semiring.
-
-Import/Export connectors to the Python ecosystem
-------------------------------------------------
-
-``graphblas.io`` contains functions for converting to and from:
-
-.. code-block:: python
-
- import graphblas as gb
-
- # numpy arrays
- # 1-D array becomes Vector, 2-D array becomes Matrix
- A = gb.io.from_numpy(m)
- m = gb.io.to_numpy(A)
-
- # scipy.sparse matrices
- A = gb.io.from_scipy_sparse(m)
- m = gb.io.to_scipy_sparse(m, format='csr')
-
- # networkx graphs
- A = gb.io.from_networkx(g)
- g = gb.io.to_networkx(A)
-
-Indices and tables
-==================
-
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
+ Find out how you can contribute and make python-graphblas
+ better for everyone.
diff --git a/docs/user_guide/collections.rst b/docs/user_guide/collections.rst
new file mode 100644
index 000000000..35b1af4a4
--- /dev/null
+++ b/docs/user_guide/collections.rst
@@ -0,0 +1,177 @@
+
+Collections
+===========
+
+python-graphblas has three collections for holding data: **Matrix**, **Vector**, and **Scalar**.
+Each is a sparse container, meaning elements may be empty. Each collection has a single data type
+and a defined shape based on the number of dimensions.
+
+Matrix
+------
+
+A **Matrix** is a 2-dimensional sparse array.
+
+.. code-block:: python
+
+ M = gb.Matrix.from_values(
+ [3, 0, 3, 5, 6, 0, 6, 1, 6, 2, 4, 1],
+ [0, 1, 2, 2, 2, 3, 3, 4, 4, 5, 5, 6],
+ [3, 2, 3, 1, 5, 3, 7, 8, 3, 1, 7, 4],
+ nrows=7, ncols=7
+ )
+
+Viewing the ``__repr__`` for ``M`` in a notebook will look like the following, although
+the name, dtype, and format might be different than what is shown here.
+
+.. image:: ../_static/img/repr-matrix.png
+
+Here are common attributes that can be accessed on the Matrix object.
+
+.. code-block:: python
+
+ >>> M.dtype
+ INT64
+ >>> M.nrows
+ 7
+ >>> M.ncols
+ 7
+ >>> M.shape
+ (7, 7)
+ >>> M.nvals
+ 12
+ >>> M.ndim
+ 2
+
+Vector
+------
+
+A **Vector** is a 1-dimensional sparse array.
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 3, 4, 6], [12.1, -5.4, 2.9, 2.2], size=8)
+
+Viewing the ``__repr__`` for ``v`` in a notebook will look like the following, although
+the name, dtype, and format might be different than what is shown here.
+
+.. image:: ../_static/img/repr-vector.png
+
+Here are common attributes that can be accessed on the Vector object.
+
+.. code-block:: python
+
+ >>> v.dtype
+ FP64
+ >>> v.size
+ 8
+ >>> v.shape
+ (8,)
+ >>> v.nvals
+ 4
+ >>> v.ndim
+ 1
+
+Scalar
+------
+
+A **Scalar** is a 0-dimensional sparse array. It either contains a single value
+or no value.
+
+.. code-block:: python
+
+ s = gb.Scalar.from_value(5)
+
+Viewing the ``__repr__`` for ``s`` in a notebook will look like the following.
+
+.. image:: ../_static/img/repr-scalar.png
+
+Here are common attributes that can be accessed on the Scalar object.
+
+.. code-block:: python
+
+ >>> s.dtype
+ INT64
+ >>> s.value
+ 5
+ >>> s.is_empty
+ False
+ >>> s.nvals
+ 1
+ >>> s.shape
+ ()
+ >>> s.ndim
+ 0
+
+Opaque Objects
+--------------
+
+Collections in GraphBLAS are opaque objects per the spec. In theory, this means that the user should not
+care how a collection is stored, leaving all those details to the backend implementation.
+
+The backend for python-graphblas is SuiteSparse:GraphBLAS, and it chooses to expose its storage formats.
+For example, a Matrix might be stored as 'csr' or 'csc' or 'bitmapr' or one of several other formats.
+
+These formats can be used to fine-tune algorithms in certain cases, but for the most part are safe to ignore.
+SuiteSparse:GraphBLAS has very good defaults which work for most use cases without any manual tuning.
+
+When manual tuning is required, look in the :ref:`api_reference` for details about the ``.ss`` namespace.
+
+Collection Methods
+------------------
+
+dup
+~~~
+
+Collections can be duplicated using ``.dup()``. This creates a fully independent object with
+the same shape, dtype, and values as the original.
+
+.. code-block:: python
+
+ w = v.dup()
+
+clear
+~~~~~
+
+Collections can be cleared using ``.clear()``. This can be used to free memory when the contents
+of an object are no longer needed, but the object will be reused later (perhaps within a loop).
+The shape and dtype remain unchanged, but the collection will be fully sparse (i.e. nvals == 0).
+
+.. code-block:: python
+
+ >>> M.clear()
+ >>> M.nvals
+ 0
+
+to_values
+~~~~~~~~~
+
+To go from a collection back to the index and values, ``.to_values()`` can be called. Numpy arrays
+will be returned in a tuple.
+
+.. code-block:: python
+
+ >>> v = gb.Vector.from_values([0, 3, 5], [2.0, 3.1, 4.7])
+ >>> idx, vals = v.to_values()
+ >>> idx
+ array([0, 3, 5], dtype=uint64)
+ >>> vals
+ array([2.0, 3.1, 4.7])
+
+Note that Scalar does not have ``.to_values()``. Instead, the ``.value`` attribute will extract
+the underlying value or return ``None`` if the Scalar is empty.
+
+Empty Constructors
+~~~~~~~~~~~~~~~~~~
+
+Empty collections can be created directly from the collection classes by specifying the dtype and dimensions.
+
+.. code-block:: python
+
+ # Create an empty 6x7 FP64 Matrix
+ M = gb.Matrix(float, 6, 7)
+
+ # Create an empty 100-element Vector with INT32 dtype
+ v = gb.Vector('int32', 100)
+
+ # Create an empty boolean Scalar
+ s = gb.Scalar(bool)
diff --git a/docs/user_guide/fundamentals.rst b/docs/user_guide/fundamentals.rst
new file mode 100644
index 000000000..f47296c47
--- /dev/null
+++ b/docs/user_guide/fundamentals.rst
@@ -0,0 +1,198 @@
+
+Fundamentals
+============
+
+The GraphBLAS specification is a mathematically sound, well thought out API. Most GraphBLAS functions
+follow a general pattern as follows:
+
+.. image:: ../_static/img/GraphBLAS-API-example.png
+
+The name of the function is obviously unique for each function, but in general the output, mask, accumulator,
+operator, input arguments, and descriptor are common to every GraphBLAS function.
+
+The **output** is the object which will be updated with the result of the computation. It must already exist
+and be initialized, even if initialized as an empty object.
+
+The **mask** is optional and, if defined, controls which elements of the output are updated. Elements outside
+the mask are left untouched.
+
+The **accumulator** is optional and, if defined, determines how existing values in the output are combined
+with new values from the computation.
+
+The **operator** is the secret sauce in GraphBLAS, allowing matrix multiplication to use a different semiring
+than the standard *plus_times*, or to add two vectors together using something other than standard addition.
+
+The **inputs** are the objects involved in the computation. There may be two inputs (for addition or matrix
+multiplication) or there may be one input (for apply or reduction).
+
+The descriptor is a set of bitwise flags.
+
+ - The first input can be transposed during the computation (only applies for Matrix input).
+ - The second input can be transposed during the computation (only applies for Matrix input).
+ - The mask can be complemented, which is needed to avoid situations where only a few elements are
+ *not* part of the mask. Rather than creating a mostly dense object, the very sparse mask can be provided
+ with a complement indicator.
+ - The structure of the mask can be used as the mask, ignoring the values in the mask object. This is more
+ efficient than checking each value and essentially treats the mask as if every element were True.
+ - Replace mode indicates that elements outside the mask area should be cleared in the final output.
+ When not in replace mode, elements outside the mask are left untouched.
+
+For more details, look at the official API spec at `graphblas.org `_.
+
+C-to-Python Mapping
+-------------------
+
+The mapping between the C API and the python-graphblas API is as follows.
+
+.. image:: ../_static/img/GraphBLAS-mapping.png
+
+All arguments which impact the output are grouped on the left side of the ``<<``.
+This includes the output itself as well as the mask, accumulator, and relevant pieces of
+the descriptor.
+
+Everything which describes the actual computation (operation, input arguments, operator)
+is grouped on the right side of the ``<<``.
+
+This separation is very intentional. While the whole expression is passed as a single call
+to the backend GraphBLAS implementation, it really is two operations fused together for
+efficiency. The computation on the right could be computed fully and subsequently the output could be
+updated with the computed object. This is obviously inefficient when a mask is used and
+many of the elements computed on the right would be thrown away. This is the reason to fuse
+the two operations for performance. However, keeping them separated in the python-graphblas
+notation makes it much easier to see and reason about the two aspects of the fused operation.
+
+The descriptor bit flags are common in C, but with Python we can have a much nicer syntax.
+The bit flags become attributes on the objects they impact.
+
+ - For input arguments, adding ``.T`` will indicate that the argument should be transposed in the operation.
+ - For the mask, ``.S`` indicates the structure is only considered, while ``.V`` indicates that the value
+ of each element should be used.
+ - The mask also has a prefix ``~`` to indicate the complement of the mask.
+ - Replace mode is indicated with a boolean keyword argument named ``replace``.
+
+Update Notation (<<)
+--------------------
+
+With the goal to separate output-affecting items on the left side of an assignment syntax, the obvious choice
+would have been to use ``M = ...``. Unfortunately, in Python this replaces what the variable ``M`` is pointing
+to rather than updating ``M``.
+
+.. code-block:: python
+
+ # This doesn't actually update M
+ M = semiring.min_plus(A @ B)
+
+
+Python does have assignment notation using square brackets. While this would work to update ``M``,
+the square bracket notation lacks the ability for keyword-arguments. Doing something like this is impossible:
+
+.. code-block:: python
+
+ # This will raise a SyntaxError
+ M[mask, accum=binary.min] = semiring.min_plus(A @ B)
+
+To have the full power of keyword arguments, the choice was made to use call syntax (i.e. parentheses) on
+the output object. Python does not allow assignment into the result of a function call. However, the
+left-shift (<<) operator can work with the result of such a call.
+
+Thus the choice was made to (ab)use the left-shift operator to indicate updates. It satisfies the visual
+feeling of computation from the right side flowing to the left. The result of the left-shift operator is
+ignored, while the object on the left is updated as a side-effect.
+
+.. code-block:: python
+
+ # This is the preferred syntax for python-graphblas
+ M(mask, accum=binary.min) << semiring.min_plus(A @ B)
+
+For anyone who objects to using the left-shift operator this way, an update method exists to perform the
+same job.
+
+.. code-block:: python
+
+ # This is equivalent to the code above
+ M(mask, accum=binary.min).update(semiring.min_plus(A @ B))
+
+Delayed Objects
+---------------
+
+In python-graphblas, most computations result in delayed objects. It is only when both sides are passed in the
+update (``<<``) that everything is brought together into a single call to the GraphBLAS backend. This gives
+the benefit of nicer Python syntax with very efficient fused GraphBLAS calls for performance.
+
+The user needs to be aware of these delayed objects when inspecting them. The ``__repr__`` is very similar,
+but will usually show something like ``MatrixExpression`` to indicate a delayed object. For small objects,
+the ``__repr__`` will usually perform the computation to show the user what the result will be. For large objects,
+this is not done to avoid an expensive computation during object inspection.
+
+Delayed objects can also be explicitly converted to real objects using the ``.new()`` method.
+
+.. code-block:: python
+
+ # This is a delayed object
+ delayed = semiring.min_plus(A @ B)
+ # This is now a real object
+ real_object = delayed.new()
+
+ # Or make it real in one go
+ M = semiring.min_plus(A @ B).new()
+
+In general, delayed objects can be passed as inputs into subsequent computations without problems.
+
+.. code-block:: python
+
+ # C is a delayed object
+ C = semiring.min_plus(A @ B)
+
+ # Passing C is okay here; D is still a delayed object
+ D = C.reduce_rowwise(monoid.plus)
+
+
+However, a delayed object cannot be used as the output object without first converting it into a real object.
+
+Functional vs. Method Call Syntax
+---------------------------------
+
+Some operations (ex. reduce) only have a method form for calling.
+
+.. code-block:: python
+
+ v << A.reduce_columnwise(monoid.times)
+
+Other operations (ex. assignment) only have a syntax form.
+
+.. code-block:: python
+
+ # Fill the 2nd column of B with all 1's
+ B[:, 1] << 1
+
+However, most operations have both a functional and a method-based calling syntax.
+For method-based calls, the operator is an argument. For functional calls, the operator
+is the calling object.
+
+.. code-block:: python
+
+ # Method form for element-wise intersection
+ C << A.ewise_mult(B, binary.plus)
+
+ # Functional form for element-wise intersection
+ C << binary.plus(A & B)
+
+Comparing Objects
+-----------------
+
+Comparing objects can be tricky because of the need to check the shape, sparsity pattern, and values
+of both objects. Two convenience methods exist to simplify comparisons.
+
+1. isequal -- for checking exact matches
+2. isclose -- for checking approximately equal for floating point dtypes
+
+By default, the types do not need to match as long as the values are equivalent. If exact dtype matching
+is required, set the ``check_dtype`` keyword argument to True.
+
+.. code-block:: python
+
+ while True:
+ w = ... # perform some computation
+ # Exit the loop once converged
+ if v.isclose(w, rel_tol=1e-03):
+ break
diff --git a/docs/user_guide/index.rst b/docs/user_guide/index.rst
new file mode 100644
index 000000000..a39bc0062
--- /dev/null
+++ b/docs/user_guide/index.rst
@@ -0,0 +1,29 @@
+
+.. _user_guide:
+
+User Guide
+==========
+
+The user guide covers fundamental concepts by topic. Starting at the beginning
+and reading through will give a new user a solid foundation for using python-graphblas.
+
+Each topic is an overview, not an exhaustive list of every method signature. For
+that level of detail, see the :ref:`api_reference`.
+
+
+Topics
+------
+
+.. toctree::
+ :maxdepth: 2
+
+ fundamentals
+ collections
+ types
+ operators
+ operations
+ init
+ io
+ udf
+ recorder
+
diff --git a/docs/user_guide/init.rst b/docs/user_guide/init.rst
new file mode 100644
index 000000000..7e1fdfa00
--- /dev/null
+++ b/docs/user_guide/init.rst
@@ -0,0 +1,58 @@
+
+GraphBLAS Initialization
+========================
+
+GraphBLAS must be initialized before it can be used. This is done with the
+``graphblas.init`` function. The backend and the mode must be specified.
+
+.. code-block:: python
+
+ import graphblas as gb
+ # Context initialization must happen before any other imports
+ gb.init('suitesparse', blocking=False)
+
+ # Now we can import other items from graphblas
+ from graphblas import binary, semiring
+ from graphblas import Matrix, Vector, Scalar
+
+
+Supported Backends
+------------------
+
+The only supported backend is current ``suitesparse``, although the plan is to support
+additional backends in the future.
+
+GraphBLAS Modes
+---------------
+
+The GraphBLAS spec contains the idea of **BLOCKING** vs. **NONBLOCKING** modes of operation.
+
+Blocking mode requires immediate execution, while nonblocking mode allows for some aspects of
+a function call to be delayed. Error checking on the inputs, compatibility with the output, etc
+must be checked immediately, but the actual computation is allowed to be delayed. The hope is
+that some GraphBLAS implementations will find optimizations if they see a string of computations
+and realize that the output of an intermediate computation is not used and there is a way to
+fuse operations together.
+
+Thus far, this mechanism hasn't been fully realized. The only known usages are with adding or
+removing a single element at a time. In nonblocking mode, this is much more efficient as the
+element additions are queued up and the object container is fully computed at the last possible
+moment.
+
+Default Initialization
+----------------------
+
+If submodules are imported before calling `graphblas.init`, a default initialization occurs.
+The default backend is ``suitesparse`` and the default mode is ``nonblocking``.
+
+If ``graphblas.init`` is called after the default initialization with different parameters,
+an error will be raised.
+
+Duplicate Initialization
+------------------------
+
+If ``graphblas.init`` is called multiple times with the same arguments, no error will be
+raised. However, if the arguments do not match, an error will be raised.
+
+There is no mechanism to change the initialization arguments after graphblas has been
+initialized the first time.
diff --git a/docs/user_guide/io.rst b/docs/user_guide/io.rst
new file mode 100644
index 000000000..fe8101674
--- /dev/null
+++ b/docs/user_guide/io.rst
@@ -0,0 +1,123 @@
+
+Input/Output
+============
+
+There are several ways to get data into and out of python-graphblas.
+
+From/To Values
+--------------
+
+Scalar has ``.from_value()`` and ``.value`` to get a Python scalar into and out of python-graphblas.
+
+Matrix and Vector, instead, have a ``.from_values()`` and a ``.to_values()`` method.
+
+``.from_values()`` takes index(es) and values as either:
+
+ - Python lists
+ - Numpy arrays
+
+If no dtype is provided, the data type is inferred from the values.
+
+If no size (Vector) or nrows/ncols (Matrix) are provided, the shape is inferred based on the largest
+index found for each dimension.
+
+``.to_values()`` returns numpy arrays. Indexes are always returned as ``uint64``, while the values
+array will match the collection dtype.
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([1, 3, 6], [2, 3, 4], float, size=10)
+
+.. csv-table::
+ :header: 0,1,2,3,4,5,6,7,8,9,10
+
+ ,2.0,,3.0,,,4.0,,,
+
+.. code-block:: python
+
+ >>> idx, vals = v.to_values()
+ >>> idx
+ array([0, 1, 3], dtype=uint64)
+ >>> vals
+ array([4., 5., 6.])
+
+
+NetworkX
+--------
+
+A python-graphblas Matrix can be created from a NetworkX graph using ``gb.io.from_networkx()``.
+
+The dtype and weight label can be specified. A list of nodes can also be provided to only pull out a subset
+of the NetworkX graph.
+
+``gb.io.to_networkx()`` goes the other direction. It always returns a DiGraph with weights labelled as "weight".
+
+
+Scipy.Sparse
+------------
+
+A python-graphblas Matrix can be created from a 2-D scipy.sparse array or matrix using
+``gb.io.from_scipy_sparse()``.
+
+``gb.io.to_scipy_sparse()`` will output a 2-D scipy.sparse array given a python-graphblas Matrix.
+The scipy.sparse format can be specified. It defaults to "csr".
+
+*Note about zero-weighted edges in scipy.sparse:* scipy.sparse makes the assumption that missing values
+and zero-weighted edges can be treated identically. There are conversions within the scipy.sparse codebase
+that drop zero-weighted edges. The conversion from python-graphblas to scipy.sparse attempts to preserve
+zero-weighted edges, but the user should be aware of the potential for errors occurring when zero-weighted
+edges are handled by scipy.sparse.
+
+
+Numpy (Dense)
+-------------
+
+While not useful for very large graphs, converting to and from small dense numpy arrays can be useful.
+
+``gb.io.from_numpy()`` will convert a 1-D array into a Vector and a 2-D array into a Matrix. When converting
+from numpy, zeros are treated as missing values.
+
+``gb.io.to_numpy()`` will convert a Vector or Matrix into the dense equivalent in numpy, filling missing
+values with zero.
+
+
+SuiteSparse Export/Import
+-------------------------
+
+For the fastest possible access to the underlying data structures in python-graphblas, SuiteSparse
+extensions exist in the ``.ss`` namespace on Matrix and Vector collections.
+
+``ss.export()`` gives options to export in a specific format or to leave the underlying data structures
+as-is for the fastest access. This can be very useful when serializing a graph and sending it across the
+wire.
+
+The result is a Python dict with numpy arrays and other metadata inside. See the :ref:`api_reference`
+for more details.
+
+.. code-block:: python
+
+ >>> A.ss.export('csr')
+ {'indptr': array([0, 2, 4, 6], dtype=uint64),
+ 'col_indices': array([0, 2, 1, 2, 0, 2], dtype=uint64),
+ 'sorted_cols': True,
+ 'nrows': 3,
+ 'ncols': 3,
+ 'is_iso': False,
+ 'format': 'csr',
+ 'values': array([10, 20, 30, 40, 50, 60])}
+
+Importing requires choosing the collection and calling the ``.ss.import_*`` methods.
+A generic ``import_any()`` method will inspect the format and call the appropriate importer.
+Otherwise, the format-specific import method name may be directly called.
+
+.. code-block:: python
+
+ >>> d = A.ss.export('csr')
+ >>> d["values"][0] = -100 # modify the serialized data
+ >>> M = gb.Matrix.ss.import_csr(**d)
+
+Note that A is unchanged in the above example.
+
+The SuiteSparse export has a ``give_ownership`` option. This performs a zero-copy
+move operation and invalidates the original python-graphblas object. When extreme speed is needed or memory is
+too limited to make a copy, this option may be needed.
diff --git a/docs/user_guide/operations.rst b/docs/user_guide/operations.rst
new file mode 100644
index 000000000..8e0f87e69
--- /dev/null
+++ b/docs/user_guide/operations.rst
@@ -0,0 +1,744 @@
+
+Operations
+==========
+
+Matrix Multiply
+---------------
+
+The GraphBLAS spec contains three methods for matrix multiplication, depending on whether
+the inputs are Matrix or Vector.
+
+ - **mxm** -- Matrix-Matrix multplications
+ - **mxv** -- Matrix-Vector multiplication
+ - **vxm** -- Vector-Matrix multiplication
+
+These three methods exist on python-graphblas collections, but the preferred approach is using
+the ``@`` symbol, which indicates matrix multiplication in Python.
+
+The default semiring for matrix multiplication is ``plus_times``, but any semiring may be used
+instead.
+
+Vectors do not have a natural row or column orientation. However, when multiplied on the left side
+of a Matrix, a Vector is treated as a 1xn row matrix. When multiplied on the right side of a Matrix,
+a Vector is treated as an nx1 column matrix.
+
+**Matrix-Matrix** Multiply Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2], [1, 2, 2, 3, 3],
+ [2., 5., 1.5, 4.25, 0.5], nrows=4, ncols=4)
+ B = gb.Matrix.from_values([0, 0, 1, 1, 2, 2, 3, 3], [1, 2, 0, 1, 1, 2, 0, 1],
+ [3., 2., 9., 6., 3., 1., 0., 5.])
+
+ # These are equivalent
+ C << A.mxm(B, op='min_plus') # method style
+ C << gb.semiring.min_plus(A @ B) # functional style
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2,3
+
+ **0**,,2.0,5.0,
+ **1**,,,1.5,4.25
+ **2**,,,,0.5
+ **3**,,,,
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,3.0,2.0
+ **1**,9.0,6.0,
+ **2**,,3.0,1.0
+ **3**,0.0,5.0,
+
+.. csv-table:: C << min_plus(A @ B)
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,11.0,8.0,6.0
+ **1**,4.25,4.5,2.5
+ **2**,0.5,5.0,
+ **3**,,,
+
+**Matrix-Vector** Multiply Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2], [1, 2, 2, 3, 3],
+ [2., 5., 1.5, 4.25, 0.5], nrows=4, ncols=4)
+ v = gb.Vector.from_values([0, 1, 3], [10., 20., 40.])
+
+ # These are equivalent
+ w << A.mxv(v, op='plus_times') # method style
+ w << gb.semiring.plus_times(A @ v) # functional style
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2,3
+
+ **0**,,2.0,5.0,
+ **1**,,,1.5,4.25
+ **2**,,,,0.5
+ **3**,,,,
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3
+
+ 10.0,20.0,,40.0
+
+.. csv-table:: w << plus_times(A @ v)
+ :class: inline
+ :header: 0,1,2,3
+
+ 40.0,170.0,20.0,
+
+**Vector-Matrix** Multiply Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3], [10., 20., 40.])
+ B = gb.Matrix.from_values([0, 0, 1, 1, 2, 2, 3, 3], [1, 2, 0, 1, 1, 2, 0, 1],
+ [3., 2., 9., 6., 3., 1., 0., 5.])
+
+ # These are equivalent
+ u << v.vxm(B, op='plus_plus') # method style
+ u << gb.semiring.plus_plus(v @ B) # functional style
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3
+
+ 10.0,20.0,,40.0
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,3.0,2.0
+ **1**,9.0,6.0,
+ **2**,,3.0,1.0
+ **3**,0.0,5.0,
+
+.. csv-table:: u << plus_plus(v @ B)
+ :class: inline
+ :header: 0,1,2
+
+ 69.0,84.0,12.0
+
+Element-wise Intersection
+-------------------------
+
+Two identically shaped collections can be intersected element-wise. Locations where only one of the
+two collections contains a value will be missing in the output.
+
+The GraphBLAS spec calls this operation `eWiseMult` because it has the same behavior as sparse
+multiplication when missing values are treated as zero. As a result, ``binary.times`` is the
+default operator for element-wise intersection.
+
+The method is named ``ewise_mult``, following the spec. The functional syntax uses the Python
+symbol for intersection ``&``.
+
+Example usage:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2], [1, 2, 0, 2, 1],
+ [2.0, 5.0, 1.5, 4.0, 0.5])
+ B = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 2, 0, 1, 1, 2],
+ [3., -2., 0., 6., 3., 1.])
+
+ # These are equivalent
+ C << A.ewise_mult(B, op='min') # method style
+ C << gb.binary.min(A & B) # functional style
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,1.5,,4.0
+ **2**,,0.5,
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,3.0,-2.0
+ **1**,0.0,6.0,
+ **2**,,3.0,1.0
+
+.. csv-table:: C << min(A & B)
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,-2.0
+ **1**,0.0,,
+ **2**,,0.5,
+
+Element-wise Union
+------------------
+
+Two identically shaped collections can perform a union element-wise. Locations where only one of the
+two collections contains a value will contain that value in the output. Where they overlap, the operator
+will compute the result.
+
+The GraphBLAS spec calls this operation ``eWiseAdd`` because it has the same behavior as sparse
+addition when missing values are treated as zero. As a result, ``binary.plus`` is the
+default operator for element-wise union.
+
+There are two methods in python-graphblas for element-wise union.
+
+ - ``ewise_add``
+
+ This is the official method based on the spec. It follows the spec by using a single value as-is when
+ there is no overlap. For binary operations which are associative (plus, and, min, etc),
+ ewise_add always gives the right answer. However, when the operation is not associative (minus, div, etc),
+ ewise_add can have suprising results.
+
+ - ``7 - 3 = 4``
+ - ``7 - missing = 7``
+ - ``missing - 7 = 7`` <-- *This might seem unexpected, but it is correct*
+
+ - ``ewise_union``
+
+ This is an extension provided by SuiteSparse:GraphBLAS. It adds a ``left_default`` and ``right_default``
+ parameter that specify what the missing value should be when there is only a single value.
+
+ If we set ``left_default=0`` and ``right_default=0``, then
+
+ - ``7 - 3 = 4``
+ - ``7 - missing = 7 - 0 = 7``
+ - ``missing - 7 = 0 - 7 = -7`` <-- *This gives the expected answer*
+
+
+The functional syntax uses the Python symbol for union ``|`` for both methods. To specify that ``ewise_union``
+should be used with the functional syntax, ``left_default`` and ``right_default`` keywords are used.
+
+**eWiseAdd** Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 0, 1, 1, 1], [0, 1, 2, 0, 2],
+ [9.0, 2.0, 5.0, 1.5, 4.0], nrows=3)
+ B = gb.Matrix.from_values([0, 0, 0, 2, 2, 2], [0, 1, 2, 0, 1, 2],
+ [4., 0., -2., 6., 3., 1.])
+
+ # These are equivalent
+ C << A.ewise_add(B, op='minus') # method style
+ C << gb.binary.minus(A | B) # functional style
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,9.0,2.0,5.0
+ **1**,1.5,,4.0
+ **2**,,,
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,4.0,0.0,-2.0
+ **1**,,,
+ **2**,6.0,3.0,1.0
+
+.. csv-table:: C << A.ewise_add(B, 'minus')
+ :class: inline
+ :header: ,0,1,2,
+
+ **0**,5.0,2.0,7.0
+ **1**,1.5,,4.0
+ **2**,6.0,3.0,1.0
+
+**eWiseUnion** Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 0, 1, 1, 1], [0, 1, 2, 0, 2],
+ [9.0, 2.0, 5.0, 1.5, 4.0], nrows=3)
+ B = gb.Matrix.from_values([0, 0, 0, 2, 2, 2], [0, 1, 2, 0, 1, 2],
+ [4., 0., -2., 6., 3., 1.])
+
+ # These are equivalent
+ C << A.ewise_union(B, op='minus', left_default=0, right_default=0) # method style
+ C << gb.binary.minus(A | B, left_default=0, right_default=0) # functional style
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,9.0,2.0,5.0
+ **1**,1.5,,4.0
+ **2**,,,
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,4.0,0.0,-2.0
+ **1**,,,
+ **2**,6.0,3.0,1.0
+
+.. csv-table:: C << A.ewise_union(B, 'minus', 0, 0)
+ :class: inline
+ :header: ,0,1,2,
+
+ **0**,5.0,2.0,7.0
+ **1**,1.5,,4.0
+ **2**,-6.0,-3.0,-1.0
+
+Extract
+-------
+
+Extraction is GraphBLAS takes a subset of a Vector or Matrix based on a set of indices.
+Extraction uses normal Python slice syntax.
+
+Extraction is not a shape-preserving operation, so indexes may be remapped during the process. For
+example, extracting indices [1, 3, 5] will yield an object with indices [0, 1, 2].
+
+If the index is a list of indices or a slice, that dimension will be preserved. If the index
+is an integer, the dimension will be collapsed.
+
+ - **Matrix[**\ *list/slice*, *list/slice*\ **] -> Matrix**
+ - **Matrix[**\ *list/slice*, *int*\ **] -> Vector**
+ - **Matrix[**\ *int*, *list/slice*\ **] -> Vector**
+ - **Matrix[**\ *int*, *int*\ **] -> Scalar**
+ - **Vector[**\ *list/slice*\ **] -> Vector**
+ - **Vector[**\ *int*\ **] -> Scalar**
+
+Vector Slice Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3, 4, 6], [10., 2., 40., -5., 24.])
+
+ w << v[:4]
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3,4,5,6
+
+ 10.0,2.0,,40.0,-5.0,,24.0
+
+.. csv-table:: w << v[:4]
+ :class: inline
+ :header: 0,1,2,3
+
+ 10.0,2.0,,40.0
+
+Matrix List Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 2, 0, 1, 0, 2],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+
+ C << A[[0, 2], :]
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,1.5,4.0,
+ **2**,0.5,,-7.0
+
+.. csv-table:: C << A[[0, 2], :]
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,0.5,,-7.0
+
+Assign
+------
+
+Assignment in GraphBLAS takes a smaller collection and expands it into a larger collection based
+on defined indices. It can be thought of as the inverse of Extract. The number of indices must match the
+shape of the input collection being assigned. However, the actual index positions refer to the location
+within the output object.
+
+Assignment uses normal Python slice syntax.
+
+Smaller rank objects can be assigned if the index is an integer rather than a list or slice.
+For example, assigning a Vector into a Matrix is possible if either the row or column index is
+an integer.
+
+Assigning a Scalar is also possible for any combination of integer, list or slice index.
+
+Matrix-Matrix Assignment Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 2, 0, 1, 0, 2],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+ B = gb.Matrix.from_values([0, 0, 1, 1], [0, 1, 0, 1]
+ [-99., -98., -97., -96.])
+
+ A[::2, ::2] << B
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,1.5,4.0,
+ **2**,0.5,,-7.0
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1
+
+ **0**,-99.0,-98.0
+ **1**,-97.0,-96.0
+
+.. csv-table:: A[::2, ::2] << B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,-99.0,2.0,-98.0
+ **1**,1.5,4.0,
+ **2**,-97.0,,-96.0
+
+Matrix-Vector Assignment Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 2, 0, 1, 0, 2],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+ v = gb.Vector.from_values([2], [-99.])
+
+ A[1, :] << v
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,1.5,4.0,
+ **2**,0.5,,-7.0
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2
+
+ ,,-99.0
+
+.. csv-table:: A[1, :] << v
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,,,-99.0
+ **2**,0.5,,-7.0
+
+Vector-Scalar Assignment Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3, 4, 6], [10, 2, 40, -5, 24])
+
+ v[:4] << 99
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3,4,5,6
+
+ 10,2,,40,-5,,24
+
+.. csv-table:: v[:4] << 99
+ :class: inline
+ :header: 0,1,2,3,4,5,6
+
+ 99,99,99,99,-5,,24
+
+Apply
+-----
+
+Apply takes an operator and applies it to every non-empty element in a collection.
+The operator can be unary, index unary, or binary.
+
+For the case of binary, an additional scalar argument must be provided as either the
+left or right argument, with the other argument being provided by the collection elements.
+
+The method name is ``apply`` and the functional form simply uses the operator as a calling
+function with the collection as the argument.
+
+**Unary** Apply Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3], [10., 20., 40.])
+
+ # These are equivalent
+ w << v.apply(gb.unary.minv)
+ w << gb.unary.minv(v)
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3
+
+ 10.0,20.0,,40.0
+
+.. csv-table:: w << minv(v)
+ :class: inline
+ :header: 0,1,2,3
+
+ 0.1,0.05,,0.025
+
+**IndexUnary** Apply Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3], [10., 20., 40.])
+
+ # These are equivalent
+ w << v.apply(gb.indexunary.index)
+ w << gb.indexunary.index(v)
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3
+
+ 10.0,20.0,,40.0
+
+.. csv-table:: w << index(v)
+ :class: inline
+ :header: 0,1,2,3
+
+ 0,1,,3
+
+**Binary** Apply Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3], [10., 20., 40.])
+
+ # These are all equivalent
+ w << v.apply('minus', right=15)
+ w << gb.binary.minus(v, right=15)
+ w << v - 15
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3
+
+ 10.0,20.0,,40.0
+
+.. csv-table:: w << v.apply('minus', right=15)
+ :class: inline
+ :header: 0,1,2,3,
+
+ -5.0,5.0,,25.0
+
+Select
+------
+
+Select takes an index unary operator and applies it to every non-missing element of a collection.
+If the result is True, the element remains in the output. If the result is False, the element
+becomes missing in the result. Thus the output is a filtered version of the original collection.
+
+Upper Triangle Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 2, 0, 2, 1, 2],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+
+ # These are equivalent
+ C << A.select('triu')
+ C << gb.select.triu(A)
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,1.5,,4.0
+ **2**,,0.5,-7.0
+
+.. csv-table:: C << select.triu(A)
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,,,4.0
+ **2**,,,-7.0
+
+Select by Value Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3, 4, 6], [10., 2., 40., -5., 24.])
+
+ # These are equivalent
+ w << v.select('>=', 5)
+ w << gb.select.value(v >= 5)
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3,4,5,6
+
+ 10.0,2.0,,40.0,-5.0,,24.0
+
+.. csv-table:: w << select.value(v >= 5)
+ :class: inline
+ :header: 0,1,2,3,4,5,6
+
+ 10.0,,,40.0,,,24.0
+
+Reduce
+------
+
+Reduction reduces the number of dimensions of a collection. A Matrix can become a Vector, and a Matrix or
+Vector can be reduced to a Scalar.
+
+When reducing a Matrix to a Vector, the reduction can be done rowwise or columnwise.
+
+A monoid or aggregator is used to perform the reduction.
+
+**Matrix-to-Vector** Columnwise Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 3, 0, 1, 0, 1],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+
+ w << A.reduce_columnwise('times')
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2,3
+
+ **0**,,2.0,,5.0
+ **1**,1.5,4.0,,
+ **2**,0.5,-7.0,,
+
+.. csv-table:: w << A.reduce_columnwise('times')
+ :class: inline
+ :header: ,0,1,2,3
+
+ ,0.75,-56.0,,5.0
+
+**Matrix-to-Scalar** Example:
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 3, 0, 1, 0, 1],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+
+ s << A.reduce_scalar('max')
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2,3
+
+ **0**,,2.0,,5.0
+ **1**,1.5,4.0,,
+ **2**,0.5,-7.0,,
+
+.. csv-table:: s << A.reduce_scalar('max')
+ :class: inline
+ :header: ,,,,
+
+ 5.0
+
+**Vector-to-Scalar** Aggregator Example:
+
+.. code-block:: python
+
+ v = gb.Vector.from_values([0, 1, 3, 4, 6], [10., 2., 40., -5., 24.])
+
+ # These are equivalent
+ s << v.reduce('argmin')
+ s << gb.agg.argmin(v)
+
+.. csv-table:: v
+ :class: inline
+ :header: 0,1,2,3,4,5,6
+
+ 10.0,2.0,,40.0,-5.0,,24.0
+
+.. csv-table:: s << argmin(v)
+ :class: inline
+ :header: ,,,
+
+ 4
+
+Transpose
+---------
+
+The transpose can either a descriptor flag set on the input of a computation or the final computation
+itself.
+
+To force the transpose to be computed by itself, use it by itself as the right-hand side of a computation.
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 3, 0, 1, 0, 2],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+
+ C << A.T
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1,2,3
+
+ **0**,,2.0,,5.0
+ **1**,1.5,4.0,,
+ **2**,0.5,,-7.0,
+
+.. csv-table:: C << A.T
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,1.5,0.5
+ **1**,2.0,4.0,
+ **2**,,,-7.0
+ **3**,5.0,,
+
+Kronecker
+---------
+
+The `Kronecker product `_ of two matrices multiplies
+every element of A (m×n) by every element of B (p×q) to create a pm×qn block matrix.
+
+The Kronecker product uses a binary operator.
+
+.. code-block:: python
+
+ A = gb.Matrix.from_values([0, 0, 1], [0, 1, 0], [1., -2., 3.])
+ B = gb.Matrix.from_values([0, 0, 1, 1, 2, 2], [1, 2, 0, 1, 0, 2],
+ [2.0, 5.0, 1.5, 4.0, 0.5, -7.0])
+
+ C << A.kronecker(B, 'times')
+
+.. csv-table:: A
+ :class: inline
+ :header: ,0,1
+
+ **0**,1.0,-2.0
+ **1**,3.0,
+
+.. csv-table:: B
+ :class: inline
+ :header: ,0,1,2
+
+ **0**,,2.0,5.0
+ **1**,1.5,4.0,
+ **2**,0.5,,-7.0
+
+.. csv-table:: C << A.kronecker(B, 'times')
+ :class: inline
+ :header: ,0,1,2,3,4,5
+
+ **0**,,2.0,5.0,,-4.0,-10.0
+ **1**,1.5,4.0,,-3.0,-8.0,
+ **2**,0.5,,-7.0,-1.0,,14.0
+ **3**,,6.0,15.0,,,
+ **4**,4.5,12.0,,,,
+ **5**,1.5,,-21.0,,,
diff --git a/docs/user_guide/operators.rst b/docs/user_guide/operators.rst
new file mode 100644
index 000000000..1dc14a93c
--- /dev/null
+++ b/docs/user_guide/operators.rst
@@ -0,0 +1,466 @@
+
+Operators
+=========
+
+Operators operate on the individual elements of collections. There are various operator classes in
+GraphBLAS, each with defined mathematical properties.
+
+Unary Operators
+---------------
+
+Unary operators map one input value to another in a possibly different domain.
+
+Example usage:
+
+.. code-block:: python
+
+ # Compute the absolute value of each element in M
+ M_abs = gb.unary.abs(M).new()
+
+Common unary operators are:
+
+ - **identity** -- returns the input unchanged
+ - **abs** -- absolute value
+ - **ainv** -- additive inverse (f(x) = -x)
+ - **minv** -- multiplicative inverse (f(x) = 1/x)
+ - **lnot** -- logical not (True -> False)
+ - **bnot** -- binary not (110010 -> 001101)
+ - **one** -- result is always 1.0 (or equivalent for the dtype)
+ - **round** -- floating-point round function
+ - **floor** -- floating-point floor function
+ - **ceil** -- floating-point ceiling function
+ - **sin** -- trigonometric sine function
+ - **cos** -- trigonometric cosine function
+ - **tan** -- trigonometric tangent function
+ - **exp** -- exponential
+ - **log** -- logarithm
+
+Unary operators are located in the ``graphblas.unary`` namespace. Additional unary operators
+registered from numpy are located in ``graphblas.unary.numpy``.
+
+Binary Operators
+----------------
+
+Binary operators take two inputs and return a single value. The input domains and output domain
+do not need to match.
+
+Output accumulators are binary operators.
+
+Example usage:
+
+.. code-block:: python
+
+ # Element-wise union, taking the max for intersecting element
+ # Accumulate the result into M via addition
+ M(accum=gb.binary.plus) << gb.binary.max(A | B)
+
+Common binary operators are:
+
+ - **pair** -- result is always 1.0 (or equivalent for the dtype)
+ - **first** -- f(a, b) = a
+ - **second** -- f(a, b) = b
+ - **min** -- min(a, b)
+ - **max** -- max(a, b)
+ - **eq** -- a == b
+ - **ne** -- a != b
+ - **gt** -- a > b
+ - **lt** -- a < b
+ - **ge** -- a >= b
+ - **le** -- a <= b
+ - **plus** -- a + b
+ - **minus** -- a - b
+ - **times** -- a * b
+ - **truediv** -- a / b
+ - **fmod** -- a % b
+ - **pow** -- a ** b
+ - **atan2** -- math.atan2(a, b)
+ - **lor** -- logical or (a | b)
+ - **land** -- logical and (a & b)
+ - **lxor** -- logical xor (a ^ b)
+ - **lxnor** -- logical xnor ~(a ^ b)
+ - **bor** -- binary or
+ - **band** -- binary and
+ - **bxor** -- binary xor
+ - **bxnor** -- binary xnor
+
+Binary operators are located in the ``graphblas.binary`` namespace. Additional binary operators
+registered from numpy are located in ``graphblas.binary.numpy``.
+
+Monoids
+-------
+
+Monoids extend the concept of a binary operator to require a single domain for all inputs and
+the output. Monoids are also associative, so the order of the inputs does not matter. And finally,
+monoids have a default identity such that ``A op identity == A``.
+
+Monoids are commonly for reductions, collapsing all elements down to a single value.
+
+Example usage:
+
+.. code-block:: python
+
+ # Sum up all non-empty elements in M
+ total = M.reduce_scalar(gb.monoid.plus).value
+
+Common monoids are:
+
+ - **any** -- return either input
+ - **min** -- min(a, b)
+ - **max** -- max(a, b)
+ - **plus** -- a + b
+ - **times** -- a * b
+ - **land** -- a & b
+ - **lor** -- a | b
+ - **lxor** -- a ^ b
+ - **lxnor** -- ~(a ^ b)
+
+Monoids are located in the ``graphblas.monoid`` namespace. Additional monoids registered from
+numpy are located in ``graphblas.monoid.numpy``.
+
+Semirings
+---------
+
+Semirings are a combination of a monoid and a binary operator. The binary operator is used for the
+"multiplication" part of a dot product, while the monoid is used for the reduction.
+
+Standard matrix multiplication uses the "plus_times" semiring.
+
+Semirings are primarily used during matrix multiplication.
+
+Example usage:
+
+.. code-block:: python
+
+ C << gb.semiring.min_plus(A @ B)
+
+Common semirings are:
+
+ - **plus_times** (standard matrix multiplication)
+ - **min_plus** (used for shortest path computations)
+ - **max_plus**
+ - **min_times**
+ - **max_times**
+ - **min_max**
+ - **max_min**
+ - **min_first**
+ - **min_second**
+ - **max_first**
+ - **max_second**
+ - **plus_min**
+ - **lor_land**
+ - **land_lor**
+
+Semirings are located in the ``graphblas.semiring`` namespace. Additional semirings registered
+from numpy are located in ``graphblas.semiring.numpy``.
+
+IndexUnary Operators
+--------------------
+
+A variant of unary operators are indexunary operators. They behave exactly like unary operators,
+but the inputs are the value, the index position(s) of that value, and an thunk parameter.
+
+For example, an IndexUnary operator applied to a Matrix would be given the value, row, and column
+of each element (plus the thunk). The operator can use all of those pieces to determine an appropriate output.
+
+IndexUnary operators are used primarily in ``select`` to filter based on the index positions.
+
+Example usage:
+
+.. code-block:: python
+
+ # Select the upper triangle
+ A_upper = gb.select.triu(A).new(name="A_upper")
+
+.. image:: ../_static/img/Matrix-A-upper.png
+
+Example usage with a thunk parameter:
+
+.. code-block:: python
+
+ # Select the upper triangle, excluding the diagonal
+ A_upper = gb.select.triu(A, 1).new(name="A_strictly_upper")
+
+.. image:: ../_static/img/Matrix-A-strictly-upper.png
+
+Defined IndexUnary operators are:
+
+ **index** -- return the vector index
+ - **rowindex** -- return the matrix row index
+ - **colindex** -- return the matrix column index
+ - **diagindex** -- return the matrix diagonal index (i.e. column - row)
+ - **tril** -- lower triangle matrix (True if column >= row)
+ - **triu** -- upper triangle matrix (True if column <= row)
+ - **diag** -- matrix diagonal (True if row == column)
+ - **offdiag** -- matrix off-diagonal (True if row != column)
+ - **indexle** -- vector index <= thunk
+ - **indexgt** -- vector index > thunk
+ - **rowle** -- matrix row index <= thunk
+ - **rowgt** -- matrix row index > thunk
+ - **colle** -- matrix column index <= thunk
+ - **colgt** -- matrix column index > thunk
+ - **valueeq** -- value == thunk
+ - **valuene** -- value != thunk
+ - **valuelt** -- value < thunk
+ - **valuele** -- value <= thunk
+ - **valuegt** -- value > thunk
+ - **valuege** -- value >= thunk
+
+IndexUnary operators are located in two places.
+
+ - ``graphblas.indexunary``
+
+ All IndexUnary operators are contained here.
+ Calling the operators in the indexunary namespace will perform an ``apply`` operation.
+
+ - ``graphblas.select``
+
+ Only the IndexUnary operators which return a boolean are contained in this namespace
+ (i.e. all except rowindex, colindex, and diagindex). Calling the operators in the
+ select namespace will perform a ``select`` operation.
+
+Aggregators
+-----------
+
+Aggregators are advanced reducers. They are similar to monoids, but do not require the same input
+and output domain. They are usually efficiently constructed recipes which require several calls to
+the backend GraphBLAS implementation, but are often thought of as a single unit of computation by
+other scientific libraries.
+
+For example, ``argmax`` reduces all the elements of a Vector to a single value, but instead of returning
+the maximum value, it returns the index of the maximum value. This requires additional work beyond what
+a simple monoid can provide.
+
+Example usage:
+
+.. code-block:: python
+
+ pos_of_largest = v.reduce(gb.agg.argmax).value
+
+Common Aggregators are:
+
+ - **count** - the number of non-empty elements
+ - **argmax** - position of the largest element
+ - **argmin** - position of the smallest element
+ - **mean** - average value (i.e. sum / count)
+ - **stdp** - population standard deviation
+ - **stds** - sample standard deviation
+ - **first** - first element
+ - **last** - last element
+
+Aggregators are located in the ``graphblas.agg`` namespace.
+
+Calling an aggregator with a collection will perform a reduction to scalar operation.
+Specifying ``rowwise=True`` or ``columnwise=True`` allows performing Matrix to
+Vector reduction rather than a full reduction to scalar.
+
+Operator Type Specialization
+----------------------------
+
+When calls are made to the backend GraphBLAS implementation, all operators are typed, meaning
+``min_FP32`` and ``min_INT64`` are different operators according to the backend.
+To avoid the user needing to worry about this detail, operator classes figure out the correct
+type variant to use automatically based on the input arguments.
+
+If desired, the user may explicitly use the typed variants of operators to force a certain behavior.
+Be aware that if the collection types do not match the operator types, the collection elements
+will be type cast in the backend using C casting rules.
+
+Example usage:
+
+.. code-block:: python
+
+ # This is the normal way to compute the minimum of a Vector
+ minval = v.reduce(gb.monoid.min).value
+
+ # This will force the FP32 version of min to be used, possibly type casting the elements
+ minvalFP32 = v.reduce(gb.monoid.min['FP32']).value
+
+
+The gb.op Namespace
+-------------------
+
+As a convenience when working with operators, a single ``graphblas.op`` namespace exists which
+combines all the operators from
+
+ - graphblas.unary
+ - graphblas.binary
+ - graphblas.monoid
+ - graphblas.semiring
+
+This facilitates writing more succinct code such as:
+
+.. code-block:: python
+
+ from graphblas import op
+
+ cur_min(accum=op.min) << op.min_plus(A @ B).reduce_rowwise(op.min)
+
+In the case of name conflicts (ex. binary.min and monoid.min), only one will exist in the the
+``graphblas.op`` namespace. However, almost all functions which require a specific kind of
+operator have a mechanism to convert from an identically named operator of a different type.
+
+For example, using monoid.min as an accumulator will automatically access binary.min instead.
+This means that using the correct name from the ``op`` namespace will almost always "just work".
+
+Infix Notation
+--------------
+
+Standard Python infix notation works in python-graphblas, but may have a specific meaning for
+each symbol. Each is detailed below.
+
+The following objects will be used to demonstrate the behavior.
+
+.. csv-table:: Vector v
+ :header: 0,1,2,3,4,5
+
+ 1.0,,2.0,3.5,,9.0
+
+
+.. csv-table:: Vector w
+ :header: 0,1,2,3,4,5
+
+ 7.0,5.2,,3.0,,2.5
+
+Operating with Scalars
+~~~~~~~~~~~~~~~~~~~~~~
+
+All infix operators involving a scalar will act only on the non-empty elements of the collection.
+
+For example, ``A + 1`` will have the same number of elements as ``A`` and will not becomes
+fully dense. In other words, missing values are *not* treated as ``0 + 1 = 1``. Instead, they are
+treated as ``missing + 1 = missing``.
+
+Addition
+~~~~~~~~
+
+Addition performs an element-wise union between collections, adding overlapping elements.
+
+.. code-block:: python
+
+ v + w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ 8.0,5.2,2.0,6.5,,11.5
+
+Subtraction
+~~~~~~~~~~~
+
+Subtraction performs an element-wise union between collections, subtracting overlapping elements
+and negating any standalone elements from the right-hand object.
+
+.. code-block:: python
+
+ v - w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ -6.0,-5.2,2.0,0.5,,6.5
+
+Multiplication
+~~~~~~~~~~~~~~
+
+Multiplication performs an element-wise intersection between collections, multiplying
+overlapping elements.
+
+.. code-block:: python
+
+ v * w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ 7.0,,,10.5,,22.5
+
+Division
+~~~~~~~~
+
+True Division ( / ) performs an element-wise intersection between collections, dividing overlapping
+elements and always results in a floating-point dtype.
+
+ - ``0 / 0 = nan``
+ - ``+x / 0 = inf``
+ - ``-x / 0 = -inf``
+
+.. code-block:: python
+
+ v / w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ 0.142857,,,1.166667,,3.6
+
+Floor Division ( // ) performs an element-wise intersection between collections, performing integer
+division on overlapping elements. For floating-point inputs, the result remains floating-point,
+but all elements are whole numbers.
+
+Dividing by zero with floor division will raise a ``ZeroDivisionError``.
+
+.. code-block:: python
+
+ v // w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ 0.0,,,1.0,,3.0
+
+Modulus
+~~~~~~~
+
+Modulus performs an element-wise intersection between collections, computing the remainder
+of dividing overlapping elements.
+
+.. code-block:: python
+
+ v % w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ 1.0,,,0.5,,1.5
+
+Power
+~~~~~
+
+Power performs an element-wise intersection between collections, computing x to
+the power of y for overlapping elements.
+
+.. code-block:: python
+
+ v ** w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ 1.0,,,42.875,,243.0
+
+Comparisons
+~~~~~~~~~~~
+
+Comparisons (==, !=, >, >=, <, <=) perform an element-wise intersection between collections,
+performing the comparison for overlapping elements. The result is always boolean.
+
+**NOTE:** to compare full equality of two collections, use ``.isequal`` or ``.isclose``
+rather than ``all(A == B)``
+
+.. code-block:: python
+
+ v > w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ False,,,True,,True
+
+.. code-block:: python
+
+ v == w
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ False,,,False,,False
diff --git a/docs/user_guide/recorder.rst b/docs/user_guide/recorder.rst
new file mode 100644
index 000000000..e37be3ad6
--- /dev/null
+++ b/docs/user_guide/recorder.rst
@@ -0,0 +1,53 @@
+
+Recorder
+========
+
+While python-graphblas dispatches to GraphBLAS C API calls in the backend, many design decisions
+have been made which make it easier to write and read python-graphblas code. As a result, it may
+not always be obvious exactly how python-graphblas code relates to the backend calls.
+
+A ``Recorder`` object exists which records each of the backend calls, making it easier to see what
+calls and how many calls are being made to the backend implementation. This can be used for debugging
+as well as for easy porting of python-graphblas algorithms into C if the need arises.
+
+The following example shows the Single-Source Shortest Path algorithm in Python.
+Because there is a Python while-loop, the Recorder will be inserted inside the
+loop to avoid viewing the calls from every iteration of the loop (which is not very helpful).
+Instead, only the calls from the last iteration will be returned.
+
+.. code-block:: python
+
+ from graphblas import op, semiring, Matrix, Vector, Recorder
+
+ # Create the graph and starting vector
+ start_node = 0
+ G = Matrix.from_values(
+ [0, 0, 1, 1, 2],
+ [1, 2, 2, 3, 3],
+ [2.0, 5.0, 1.5, 4.25, 0.5],
+ nrows=4, ncols=4)
+ v = Vector.from_values([start_node], [0.0], size=4)
+
+ # Compute SSSP
+ while True:
+ # Place the recorder inside to loop to only get one-loop worth of calls
+ with Recorder() as r:
+ w = v.dup()
+ v(op.min) << semiring.min_plus(v @ G)
+ if v.isequal(w):
+ break
+
+.. image:: ../_static/img/Recorder-output.png
+
+To explain the result,
+
+ - v.dup() became a single ``GrB_Vector_dup`` call
+ - The min_plus matrix multiplication became a single ``GrB_vxm`` call
+ - isequal is a recipe and became 3 lines of backend calls
+
+ - ``GrB_Vector_new`` to create a temporary object
+ - ``GrB_Vector_eWiseMult_BinaryOp`` using the ``EQ`` operator
+ - ``GrB_Vector_reduce_BOOL`` using the ``LAND`` (logical and) operator to reduce
+
+The names of the several objects (v_70, v66, s_14) are obviously random, but with these lines,
+it should be very simple to write an equivalent SSSP algorithm using a C while-loop and C API calls.
diff --git a/docs/user_guide/types.rst b/docs/user_guide/types.rst
new file mode 100644
index 000000000..ea7ca0824
--- /dev/null
+++ b/docs/user_guide/types.rst
@@ -0,0 +1,45 @@
+
+Data Types
+==========
+
+Each collections must have a single data type, indicated by the ``.dtype`` attribute.
+
+Defined data types are:
+
+ - BOOL
+ - UINT8
+ - UINT16
+ - UINT32
+ - UINT64
+ - INT8
+ - INT16
+ - INT32
+ - INT64
+ - FP32
+ - FP64
+ - FC32 (complex float32, not supported on Windows)
+ - FC64 (complex float64, not supported on Windows)
+
+Each of these defined types has a string representation (accessed by ``.name``),
+a corresponding dtype in numpy (accessed by ``.np_type``), and a corresponding
+dtype in numba (accessed by ``.numba_type``).
+
+When API calls need a dtype, a string or numpy or numba dtype may be used.
+Additionally, the Python builtin ``bool``, ``int``, and ``float`` may be used.
+``int`` indicates INT64 and ``float`` indicates FP64.
+
+The ``graphblas.dtypes`` namespace contains objects for each of the dtypes.
+
+User-defined Types
+------------------
+
+python-graphblas supports user-defined types.
+
+First create a custom numpy dtype. Then register it in the ``graphblas.dtypes`` namespace.
+
+.. code-block:: python
+
+ NP_Point = np.dtype([("x", np.int64), ("y", np.int64)], align=True)
+ Point = gb.dtypes.register_new("Point", NP_Point)
+ # Create a 10-element sparse vector holding Points
+ v = gb.Vector(Point, size=10)
diff --git a/docs/user_guide/udf.rst b/docs/user_guide/udf.rst
new file mode 100644
index 000000000..37923b4bd
--- /dev/null
+++ b/docs/user_guide/udf.rst
@@ -0,0 +1,54 @@
+
+User-defined Functions
+======================
+
+1. Show example of register_new
+2. Discuss commonality for other operators
+3. Discuss register_anonymous
+
+Python-graphblas requires ``numba`` which enables compiling user-defined Python functions
+to native machine code for use by the GraphBLAS backend. This provides functions which are
+very performant.
+
+Example user-defined UnaryOp:
+
+.. code-block:: python
+
+ from graphblas import unary
+ from graphblas.operator import UnaryOp
+
+ def force_odd_func(x):
+ if x % 2 == 0:
+ return x + 1
+ return x
+
+ UnaryOp.register_new('force_odd', force_odd_func)
+
+ v = Vector.from_values([0, 1, 3, 4, 5], [1, 2, 3, 8, 14])
+ w = v.apply(unary.force_odd).new()
+
+.. csv-table:: w
+ :header: 0,1,2,3,4,5
+
+ 1,3,,3,9,15
+
+
+Similar methods exist for BinaryOp and IndexUnaryOp. User-defined Monoids and Semirings are
+constructed out of previously defined and built-in UnaryOps and BinaryOps.
+
+Auto-registration of Lambdas
+----------------------------
+
+As a convenience, any lambda expression used in place of a UnaryOp will be automatically
+compiled as registered anonymously.
+
+Example lambda usage:
+
+.. code-block:: python
+
+ v.apply(lambda x: x % 5 - 2).new()
+
+.. csv-table::
+ :header: 0,1,2,3,4,5
+
+ -1,0,,1,1,2
diff --git a/graphblas/operator.py b/graphblas/operator.py
index 9e04f95eb..49ecbab44 100644
--- a/graphblas/operator.py
+++ b/graphblas/operator.py
@@ -1241,7 +1241,7 @@ class IndexUnaryOp(OpBase):
"trim_from_front": 4,
"num_underscores": 1,
"re_exprs": [
- re.compile("^GrB_(ROWINDEX|COLINDEX)_(INT32|INT64)$"),
+ re.compile("^GrB_(ROWINDEX|COLINDEX|DIAGINDEX)_(INT32|INT64)$"),
],
"re_exprs_return_bool": [
re.compile("^GrB_(TRIL|TRIU|DIAG|OFFDIAG|COLLE|COLGT|ROWLE|ROWGT)$"),
diff --git a/notebooks/Intro to GraphBLAS + SSSP example.ipynb b/notebooks/Intro to GraphBLAS + SSSP example.ipynb
index 9d9de06ba..ac08d6a80 100644
--- a/notebooks/Intro to GraphBLAS + SSSP example.ipynb
+++ b/notebooks/Intro to GraphBLAS + SSSP example.ipynb
@@ -366,7 +366,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
@@ -380,7 +380,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.8.3"
+ "version": "3.9.7"
}
},
"nbformat": 4,