diff --git a/.gitignore b/.gitignore index 684d34fadc1f13650e44cea8d67d64e10015354f..3453fe68f715b46d48fb1b2758c90fb58c08f202 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,15 @@ *~ __pycache__ *.egg-info + +!doc/ +!doc/Makefile + +!doc/source/*.rst +!doc/source/_static +!doc/source/_templates + +!lappy/lib +!lappy/target + +lappy/version.py diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4c09dc79472c1218dc7814e7c505b5fd95ac84cf..bb102b21581f73f8cdb51801a7643c2893b7f45f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -27,8 +27,34 @@ Python 3 POCL: Flake8: script: - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-flake8.sh - - ". ./prepare-and-run-flake8.sh lappy test" + - ". ./prepare-and-run-flake8.sh lappy test experiments" tags: - python3 except: - tags + +Documentation Sphinx: + script: + - export PY_EXE=python3 + - export PYOPENCL_TEST=portable + - export EXTRA_INSTALL="pybind11 numpy sphinx sphinx_rtd_theme" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project.sh + - ". ./build-py-project.sh" + - cd doc && make html + - git config --global user.name "${COMMIT_USER}" + - git config --global user.email "${COMMIT_EMAIL}" + - git clone https://${GH_USER}:${GH_ACCESS_TOKEN}@${HOMEPAGE_URL} homepage + - cd homepage + - git checkout src + - cd static && mkdir -p docs + - rm -rf docs/lappy + - cp -r ../../build/html docs/lappy + - git add -f ./docs/lappy + - (git commit -m "Auto-updated lappy docs") || (echo "Docs are up to date") + - git push + tags: + - python3 + - pocl + except: + - tags + allow_failure: true diff --git a/README.rst b/README.rst index e8c921c34d878111e5825002d463b5f3e9a3ad91..e3314d3d0ad429c9cd8f073a9c92861140a8c7bb 100644 --- a/README.rst +++ b/README.rst @@ -8,6 +8,11 @@ An array object has six basic components: the shape, the dtype, the expression, Scalars are special arrays with shape being an empty tuple. If the dtype is an integer type, there can be assumptions associated with the array elements represented by an integer set (for scalars) or map (for non-scalars). +Documentation +------------- + +Please visit `https://xiaoyu-wei.com/docs/lappy/`_ for the latest documentation. + Status ------ diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..d0c3cbf1020d5c292abdedf27627c6abe25e2293 --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/doc/make.bat b/doc/make.bat new file mode 100644 index 0000000000000000000000000000000000000000..6247f7e231716482115f34084ac61030743e0715 --- /dev/null +++ b/doc/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/doc/source/_static/.gitkeep b/doc/source/_static/.gitkeep new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/doc/source/_templates/.gitkeep b/doc/source/_templates/.gitkeep new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/doc/source/conf.py b/doc/source/conf.py new file mode 100644 index 0000000000000000000000000000000000000000..ce73d9ef20453c83a4d245fad06a1f9fc00a38ec --- /dev/null +++ b/doc/source/conf.py @@ -0,0 +1,72 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +sys.path.insert(0, os.path.abspath('..')) + +import lappy + +# -- Project information ----------------------------------------------------- + +project = 'Lappy' +copyright = lappy.__copyright__ +author = 'Lappy Developers' + +# The full version, including alpha/beta/rc tags +release = lappy.__version__ +version = lappy.version.GIT_REVISION + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.todo', + 'sphinx.ext.autosummary', + 'sphinx.ext.coverage', + 'sphinx.ext.mathjax', + 'sphinx.ext.ifconfig', + 'sphinx.ext.viewcode', + 'sphinx.ext.intersphinx', + 'sphinx.ext.githubpages'] + +templates_path = ['_templates'] +master_doc = 'index' +source_suffix = '.rst' + +mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_CHTML" # noqa: E501 + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = True + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] diff --git a/doc/source/index.rst b/doc/source/index.rst new file mode 100644 index 0000000000000000000000000000000000000000..4c7378c11e6c80d46a7dc6b2694650e8fd73f83b --- /dev/null +++ b/doc/source/index.rst @@ -0,0 +1,23 @@ +.. Lappy documentation master file, created by + sphinx-quickstart on Tue Apr 7 09:42:33 2020. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to Lappy's documentation! +================================= + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + ndarray + loop_domain + + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/doc/source/loop_domain.rst b/doc/source/loop_domain.rst new file mode 100644 index 0000000000000000000000000000000000000000..7943308f8cb7fabb5450c17711b22880c4df17f7 --- /dev/null +++ b/doc/source/loop_domain.rst @@ -0,0 +1,214 @@ +Loop Domain +=========== + +.. note:: + + (Vaporware Warning). + This section mostly contain pseudo-codes that are not implemented yet. + The plans may change at anytime. + +The loop domain of a lazy array is an integer set of all work items. +In :mod:`lappy`, it is represented as a set expression. + +Expression Tree +--------------- + +The expression tree of an array ``domain_expr`` starts from leaf nodes +of type :class:`SetVar`, :class:`SetParam` or :class:`SetZero`. +It has two levels: + +1. (Affine Expressions). + The nodes close to the leaf level that are of type :class:`PwAff`. + They consist the expression trees of the (irreducible) + low-dimensional polygons of type :class:`PwAffComparison`. + To avoid branches, the polygons' subtrees should have disjoint + variables (inames), but may share parameters (loop bounds). + +2. (Integer Set Expression). + The remaining nodes of type :class:`SetIntersection` and :class:`SetUnion` + represent the expression from the basic polygons + to the loop domain. + +The integer set expression exposes structures within the loop domain +that could benefit array manipulations, e.g., + +1. The :class:`PwAff` objects are labeled consistently. ISL treat inames + as purely positional placeholders, which may not be as handy. For example, + an FEM code may want to tag an iname as ``'i_q_point'``. Position of + the iname may change as the array is manipulated (e.g. transpose). + +2. The affine expressions are nonlinearity-tolerant w.r.t paramters. + ISL only allow affine operations among paramters. + When assembling the array expression, some loop bounds do not become + affine until more shaped information are fixed (e.g. reshape). + + +Temporaries and CSEs +-------------------- + +Here is a sample integer set expression, where +``OR`` represents set union, ``AND`` represents set intersection. + +Consider the following pseudo-code: + +.. code-block:: python + + n = Unsigned('n') + k = arange(n) + c = outer(sin(k), cos(k)) + +The expression for ``c`` is flattened by default, yielding a single instruction like +(names are replaced for readability). + +.. code-block:: + + for i, j + c[i, j] = sin(i) * cos(j) + end + +If the user choose to make ``k`` a temporary, + +.. code-block:: python + + n = Unsigned('n') + k = arange(n) + k.pin() + c = outer(sin(k), cos(k)) + +The instructions for ``c`` shall be like + +.. code-block:: + + for i, j + <> k[i] = i {id=k} + c[i, j] = sin(k[i]) * cos(k[j]) {dep=k} + end + +The loop domain is unchanged, +but depending on the loop nesting order, ``k`` may be computed multiple times. + +.. code-block:: + + AND + /-------/ \-------\ + | | + { 0 <= i < n } { 0 <= j < n } + +If the user choose to make ``k`` a CSE, + +.. code-block:: python + + n = Unsigned('n') + k = arange(n) + k.pin_cse() + c = outer(sin(k), cos(k)) + +The inames and temporaries of ``k`` are then duplicated. +It should produce the same results as evaluating ``k`` separately +by calling ``k.eval()`` and pass in its value as an argument. + +.. code-block:: + + for i_dup + <> k[i_dup] = i_dup {id=k} + end + + for i, j + c[i, j] = sin(k[i]) * cos(k[j]) {dep=k} + end + + +And the loop domain acquires a new branch, + +.. code-block:: + + OR + /------------------/ \-------------\ + | | + AND { 0 <= i_dup < n } + /-------/ \-------\ + | | + { 0 <= i < n } { 0 <= j < n } + + + +Diamond Problem +--------------- + +The following code produces a diamond structure in the computation graph: + +.. code-block:: python + + A = ndarray('A', shape='n, n') + x = ndarray('x', shape='n') + k = A @ x + sk = sin(k) + ck = cos(k) + C = outer(sk, ck) + +The data flow looks like (top to bottom): + +.. code-block:: + + A @ x + | + k + /\ + / \ + sk ck + \ / + \/ + c + +By default, the lazy evaluation expands all expressions, yielding +redundant computations. + +.. code-block:: + + c[i, j] = sin(sum([k], A[i, k] * x[k])) * cos(sum([k], A[j, k] * x[k])) + +Note that the ``matvec`` is computed twice above. + +Besides Loopy's ability to recognize optimization opportunities, user can +also make ``k`` a temporary or a CSE. + +.. code-block:: python + + A = ndarray('A', shape='n, n') + x = ndarray('x', shape='n') + k = A @ x + k.pin() + sk = sin(k) + ck = cos(k) + C = outer(sk, ck) + +Making ``k`` a temporary yields: + +.. code-block:: + + <> k[i] = sum([k], A[i, k] * x[k]) {id=k} + c[i, j] = sin(k[i]) * cos(k[j]) {dep=k} + +Similarly, + +.. code-block:: python + + A = ndarray('A', shape='n, n') + x = ndarray('x', shape='n') + k = A @ x + k.pin_cse() + sk = sin(k) + ck = cos(k) + C = outer(sk, ck) + +Making ``k`` a CSE yields: + +.. code-block:: + + for i_dup + <> k[i_dup] = sum([k], A[i_dup, k] * x[k]) {id=k} + end + + for i, j + c[i, j] = sin(k[i]) * cos(k[j]) {dep=k} + end diff --git a/doc/source/ndarray.rst b/doc/source/ndarray.rst new file mode 100644 index 0000000000000000000000000000000000000000..daed0d4913c534dca5a3320bea211efdc0e841c2 --- /dev/null +++ b/doc/source/ndarray.rst @@ -0,0 +1,4 @@ +NDArray +======= + +.. automodule:: lappy.ndarray diff --git a/experiments/finite_difference.py b/experiments/finite_difference.py new file mode 100644 index 0000000000000000000000000000000000000000..4406f685d0da1ef41f877e67fb4fa1b36e09e382 --- /dev/null +++ b/experiments/finite_difference.py @@ -0,0 +1,95 @@ +"""One-dimensional advection-diffusion by the FTCS scheme + + df/dt + U df/fx = D d2f/dx2 +""" +import numpy as np +import pyopencl as cl + +import lappy as lap +import lappy.core.ufuncs +import lappy.core.array + +from lappy import to_lappy_array + + +def solve_analytic( + nx=200, nt=2, dt=1e-4, + U=1., D=0.05, length=2.): # noqa + """Exact solution. + """ + final_time = dt * nt + h = length / (nx - 1) + xx = np.arange(nx) * h + + sin = np.sin + exp = np.exp + + def exact(time): + ex = exp(-4 * np.pi * np.pi * D * time) * ( + 0.5 * sin(2 * np.pi * (xx - time))) + return ex + + return exact(final_time) + + +def solve_numpy( + nx=200, nt=2, dt=1e-4, + U=1., D=0.05, length=2.): # noqa + + h = length / (nx - 1) + xx = np.arange(nx) * h + + # initial conditions + f = 0.5 * np.sin(2 * np.pi * xx) + + for m in range(nt): + f0 = f + fp = np.roll(f0, -1) + fm = np.roll(f0, 1) + + f = ( + f0 + - U * 0.5 * (dt / h) * (fp - fm) + + D * (dt / h**2) * (fp - 2 * f0 + fm) + ) + return f + + +def solve_lappy( + nx=200, nt=2, dt=1e-4, + U=1., D=0.05, length=2.): # noqa + + h = length / (nx - 1) + xx = np.arange(nx) * h + + env = lap.core.array.default_env() + env['cl_ctx'] = cl.create_some_context() + xx = to_lappy_array(xx, base_env=env) + + def roll(x, y): + # FIXME + return x + + # initial conditions + f = 0.5 * lap.sin(2 * np.pi * xx) + + for m in range(nt): + f0 = f + fp = roll(f0, -1) + fm = roll(f0, 1) + + f = ( + f0 + - U * 0.5 * (dt / h) * (fp - fm) + + D * (dt / h**2) * (fp - 2 * f0 + fm) + ) + + return f + + +if __name__ == '__main__': + fe = solve_analytic() + fn = solve_numpy() + fl = solve_lappy() + +# flake8: noqa diff --git a/experiments/lazy_arange.py b/experiments/lazy_arange.py new file mode 100644 index 0000000000000000000000000000000000000000..2baab520a80ea2d76ea378e6c32580ea5be17e80 --- /dev/null +++ b/experiments/lazy_arange.py @@ -0,0 +1,17 @@ +import pyopencl as cl +from lappy.core.function_base import arange + +ctx = cl.create_some_context() +queue = cl.CommandQueue(ctx) + +a = arange(64) +# print(a.eval(queue)) + +b = a.reshape([2, 2, 16]) +# print(b.eval(queue)) + +# FIXME: there needs to be a barrier after reshape's temporary + +c = b[:, 1, ::2] +print(c.ndim) +print(c.eval(queue)) diff --git a/experiments/multiinsn.py b/experiments/multiinsn.py new file mode 100644 index 0000000000000000000000000000000000000000..a917c26e385fad267575a0b7abff8592c7fc452b --- /dev/null +++ b/experiments/multiinsn.py @@ -0,0 +1,13 @@ +import pyopencl as cl +from lappy.core.function_base import arange + +ctx = cl.create_some_context() +queue = cl.CommandQueue(ctx) + +a = arange(10).with_name('A') +a.pin() + +b = a + 1 +# print(a.eval(queue)) +print(b.eval(queue)) +print(a.namespace[a.name][1]) diff --git a/lappy/__init__.py b/lappy/__init__.py index f6e08c472815031240ba8d1b75ff2be040e8d18f..2ab68458402462699e51ad42a29583b15da87798 100644 --- a/lappy/__init__.py +++ b/lappy/__init__.py @@ -22,7 +22,42 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from lappy.core import Array, to_lappy_array, to_lappy_shape +__all__ = [] +from lappy.version import GIT_REVISION as __git_revision__ # noqa: N811 +from lappy.version import VERSION_TEXT as __version__ # noqa: N811 +__all__ += ['__git_revision__', '__version__'] -__all__ = ["Array", "to_lappy_array", "to_lappy_shape"] +from lappy.exceptions import AxisError, ComplexWarning, RankWarning +__all__ += ['AxisError', 'ComplexWarning', 'RankWarning'] + +from lappy.core import to_lappy_array +from lappy.core.tools import ScalarType +__all__ += ['to_lappy_array', 'ScalarType'] + +from lappy.core.ufuncs import sin, cos, exp +__all__ += ['sin', 'cos', 'exp'] + +from lappy.core.array import Array +from lappy.eval import Compiler, Executor + + +class ndarray(Array): # noqa: N801 + def eval(self, queue, check_preconditions=False): + if self.value is None: + if hasattr(self, '_closure'): + pass + else: + compiler = Compiler(check_preconditions) + self._closure = compiler(self) + + print(self._closure.kernel) + print(self._closure.data_map) + + evaluator = Executor(self._closure) + res = evaluator(queue) + self.value = res[self.name] + return self.value + + +__all__ += ["ndarray", ] diff --git a/lappy/core/__init__.py b/lappy/core/__init__.py index e3b5af4235ec26a105f84dc5ae9e5a25311bf3ca..ff6cbde5a5ee9e8e14483e96581fb45437e49c8b 100644 --- a/lappy/core/__init__.py +++ b/lappy/core/__init__.py @@ -22,7 +22,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from lappy.core.array import Array, to_lappy_array, to_lappy_shape +from lappy.core.array import ( + Array, to_lappy_array, to_lappy_shape) __all__ = ["Array", "to_lappy_array", "to_lappy_shape"] diff --git a/lappy/core/array.py b/lappy/core/array.py index 3c9e95718130d6d80e13162b3523a731e7e277b7..08b0aa284f0c075108da1a61b55b15131c7ae1b5 100644 --- a/lappy/core/array.py +++ b/lappy/core/array.py @@ -1,6 +1,17 @@ from __future__ import division, absolute_import, print_function -copyright__ = "Copyright (C) 2017 Sophia Lin and Andreas Kloeckner" +import warnings +import numpy as np + +import islpy as isl +import pymbolic +from pymbolic import var, evaluate +from pymbolic.primitives import Lookup, Subscript, Variable, Product +from pymbolic.mapper.evaluator import UnknownVariableError + +from pprint import pformat + +copyright__ = "Copyright (C) 2019 Sophia Lin, Andreas Kloeckner and Xiaoyu Wei" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -21,55 +32,30 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import warnings -import logging -from traceback import format_list, extract_stack -from itertools import chain - -import pyopencl as cl - -import numpy as np -import loopy as lp - -from operator import mul -import functools - -import islpy as isl -import pymbolic -from pymbolic import var, evaluate -from pymbolic.primitives import ( - Expression, Substitution, Lookup, Subscript, Variable, Product) -from pymbolic.mapper.evaluator import UnknownVariableError - -from pprint import pformat def default_env(): """Default environment to be captured. """ - return { - "isl_ctx": isl.DEFAULT_CONTEXT, - "cl_ctx": None, - "cu_ctx": None - } + return dict() -class PreconditionNotMetError(Exception): - pass - # {{{ base lazy obj class class LazyObjectBase(object): """Base class for objects that are lazily evaluated. - Lazy objects are identified by their names. If not otherwise - specified, a unique name is generated on construction. + Lappy's lazy objects (Lapjects) are identified by their names. + If not otherwise specified, a unique name is generated. .. attribute:: name the name of the object, should be a Python identifier and do not start with double underscore '__' (which is reserved - for internal use). + for internal use). A lazy object's symbolic information is + uniquely determined by the name, that is, if two lazy objects + share the same name, they are treated as the same computation + (but can have different captured environments). .. attribute:: expr @@ -77,76 +63,236 @@ class LazyObjectBase(object): .. attribute:: value - the actual data. None if not yet evaluated + the actual data. None if not yet evaluated. + + .. attribute:: env - .. attribute:: arguments + the captured environment, a :class:`dict`. The environment is + fully known before runtime, i.e., does not contain lazy objects. - the formal argument list, a :class:`dict`, with keys :class:`str` - of lazy object names and values the corresponding lazy objects - that are stateless and whose values are not computable. + If env is None, the object is called **stateless**. Stateless objects serve + as parts of a lazy object with state, e.g., as the input arguments. - As a special case, when the argument means "self", to avoid infinite - recursion, the value should be set to None. + .. attribute:: namespace - .. attribute:: bound_arguments + a dictionary with names as keys, and tuple of - the arguments that are bound to an input variables (e.g. a numpy - array). The data is looked up in the captured environment by its - name (bound arguments are also stateless). + ``(stateless_lazy_object, object_desc)`` - .. attribute:: intermediaries + as values, where ``object_desc`` is a dictionary of tags. Sample tags are: - (conceptual) intermediate lazy objects (objects produced in the - process of making the array expression). Those objects should be stored - as stateless and shares the same env as the current array. + ``is_argument``: whether the object is an input argument. + ``is_temporary``: whether the object is a temporary. - .. attribute:: env + /All tags default to False./ - the captured environment, a :class:`dict`. The environment must be - fully known before runtime, i.e., cannot contain lazy objects. + The namespace associated with a lazy object contains all + lazy objects involved during the construction of itself. Stateless + lazy objects' namespaces should be None. - If env and value are both None, the object is called **stateless**. + In particular, to avoid duplications, ``None`` is used wherever infinite + recursion occurs, e.g. ``self.namespace[self.name][0] = None``. + + Besides the input arguments, the namespace can also contain temporaries + produced in the process of assembling the expression. Those temporaries do + not translate to actual temporary variable in the kernel. But they could + contain useful information otherwise lost, e.g. user-specified names of the + intermediate results. .. attribute:: preconditions preconditions are checked at runtime, prior to calling the kernel. It is a list of callables, each of which can take one or more - arguments/bound_arguments. If a precondition fails, an - :class:`Exception` is thrown (by either the checker - function itself, or by returning False from the checker). + input arguments. If a precondition fails, an + :class:`Exception` is thrown prior to calling the kernel. """ _counter = 0 _name_prefix = '__lappy_object_' + # {{{ constructor + def __init__(self, **kwargs): - """The constructor is not intended for user invokation.""" + """Constructor. + + :param name: (optional) name of the lapject + :param expr: the expression + :param env: (optional) captured environment + :param namespace: (optional) list of symbols + :param value: (optional) value of the lapject + :param preconditions: dynamic pre-runtime checks + """ + self.name = kwargs.pop("name", None) self.expr = kwargs.pop("expr") - self.bound_arguments = kwargs.pop("bound_arguments", dict()) - self.intermediaries = kwargs.pop("intermediaries", dict()) - self.env = kwargs.pop("env", default_env()) - self.preconditions = kwargs.pop("preconditions", list()) if self.name is None: self.name = self._make_default_name() - # only set value when given (which overwrites the value given in env) + self.namespace = kwargs.pop("namespace", None) + self.env = kwargs.pop("env", None) + + if self.namespace is None: + self.namespace = {self.name: (None, {'is_argument': True})} + + if self.name not in self.namespace: + self.namespace[self.name] = (None, {'is_argument': True}) + if 'value' in kwargs.keys(): - self.value = kwargs.pop("value") + if self.is_stateless: + raise TypeError("cannot construct a stateless obj with value") + elif self.value is not None: + if self.value != kwargs['value']: + raise ValueError( + "got conflicting values set from env and value") + else: + pass + else: + self.value = kwargs.pop("value") + + self.preconditions = kwargs.pop("preconditions", list()) + + # }}} End constructor + + # {{{ namespace/env utils + + def _decl(self, lapject, tags=None): + """Add symbol to the namespace. + """ + if self.is_stateless: + raise TypeError( + "can only declare names inside a stateful " + "lappy object ") + assert isinstance(lapject, LazyObjectBase) + name = lapject.name + + if name in self.namespace: + raise ValueError("name %s is taken" % name) + + if lapject.has_state: + self._meld_namespace(lapject.namespace, lapject.as_stateless()) + self._meld_env(lapject.env) + if tags is None: + tags = lapject.namespace[name][1].copy() + + if tags is None: + tags = dict() + else: + assert isinstance(tags, dict) + + self.namespace[name] = (lapject.as_stateless(), tags) + + def _check_namespace_integrity(self, namespace=None, name=None): + """Check that all namespace keys equals the name of the contained + object. + """ + if namespace is None: + namespace = self.namespace + if name is None: + name = self.name + for key, (val, tag) in namespace.items(): + if val is None: + assert name == key + else: + assert val.name == key - self.arguments = kwargs.pop("arguments", {self.name: None}) + def _meld_namespace(self, namespace, stateless_self=None): + """Merge with another namespace. The other namespace may contain a + self-reference (reference to the name containing the namespace). + If so, the ``stateless_self`` is used to substitute ``None``. + + If there are multiple definitions of a name, the function will + try to merge the information from all sources. An exception is + raised if conflicting information is found. + """ + handled_selref = False + for nm, (lapject, tags) in namespace.items(): + if lapject is None: + if stateless_self is None: + raise ValueError("a stateless copy of the enclosing " + "lapject is needed when importing " + "a namespace with self-reference") + if handled_selref: + raise ValueError("illegal namespace (contains more " + "than one self-references)") + assert isinstance(stateless_self, LazyObjectBase) + assert stateless_self.is_stateless + obj = stateless_self + handled_selref = True + else: + obj = lapject + + if nm in self.namespace: # try to merge information + sfobj = self.namespace[nm][0] + sftags = self.namespace[nm][1] + if sfobj is None: + sfobj = self + if sfobj.ndim != obj.ndim: + raise ValueError("conflicting ndims") + if sfobj.shape != obj.shape: + raise ValueError("conflicting shapes") + if sfobj.inames != obj.inames: + raise ValueError("conflicting inames") + + # merge tags + for tag, tval in tags.items(): + if tag in self.namespace[nm][1]: + if tval != sftags[tag]: + raise ValueError( + "conflicting tag values: %s" % tag) + else: + sftags[tag] = tval + + else: + self.namespace[nm] = (obj.as_stateless(), tags) + + def _meld_env(self, env): + """Merge with another env. + """ + for key, val in env.items(): + if key in self.env: + if self.env[key] is val: + pass # same value + elif val is None: + pass + elif self.env[key] is None: + self.env[key] = val + else: + raise ValueError( + "Trying to merge environments with conflicting " + "definitions of %s" % key) + else: + self.env[key] = val + + # }}} End namespace/env utils + + def _extract(self, member_name): + """Given a name, returns the lapject from the namespace, with the state + attached. + """ + assert isinstance(member_name, str) + if member_name not in self.namespace: + raise ValueError("name %s not found" % member_name) + if member_name == self.name: + return self + lapject = self.namespace[member_name][0].with_state( + env=self.env.copy(), namespace=self.namespace.copy()) + lapject.namespace[self.name] = ( + self.as_stateless(), + self.namespace[self.name][1].copy()) + lapject.namespace[member_name] = ( + None, + self.namespace[member_name][1].copy()) + return lapject def _make_default_name(self): - name = '%s%d' % (self.__class__._name_prefix, self.__class__._counter) - self.__class__._counter += 1 - return name + return make_default_name(self.__class__) def __str__(self): """Displays the value if available, otherwise display the type and name for user-defined names or the expression for the auto-generated names. """ - if self.value is not None: + if (not self.is_stateless) and (self.value is not None): return str(self.value) if len(self.name) >= 2 and self.name[:2] == '__': # auto-generated name @@ -168,7 +314,8 @@ class LazyObjectBase(object): elif isinstance(val, tuple): repr_dict[key] = '(' + ', '.join(str(v) for v in val) + ')' elif isinstance(val, dict): - repr_dict[key] = '{' + ', '.join(str(v) for v in val.keys()) + '}' + repr_dict[key] = '{' + ', '.join( + str(v) for v in val.keys()) + '}' else: repr_dict[key] = str(val) @@ -188,34 +335,91 @@ class LazyObjectBase(object): def __eq__(self, other): """Two objects are considered equal if they share the same name. - Also allows comparing with scalar constants, which amounts to comparing the - value with other. + Also allows comparing with scalar constants, which amounts to + comparing the value with other. """ if np.isscalar(other): return self.value == other else: return self.name == other.name + def __bool__(self): + """If then value is not known, the boolean value is always false. + (Think of it as a simplified trinary logic). + """ + if self.value is None: + return False + + return bool(self.value) + + @property + def is_stateless(self): + return self.env is None + + @property + def has_state(self): + return self.env is not None + + @property + def arguments(self): + """The input argument list. + """ + if self.is_stateless: + raise TypeError("stateless lapject has no arguments") + return { + key: val[0] + for key, val in self.namespace.items() + if val[1].get('is_argument', True) and ( + key not in self.env) + } + + @property + def bound_arguments(self): + """The arguments that are bound to data. + """ + if self.is_stateless: + raise TypeError("stateless lapject has no arguments") + return { + key: val[0] + for key, val in self.namespace.items() + if val[1].get('is_argument', True) and ( + key in self.env) + } + + @property + def temporaries(self): + """The temporaries (internal nodes of the computation graph). + """ + if self.is_stateless: + raise TypeError("stateless lapject has no temporaries") + + return { + key: val[0] + for key, val in self.namespace.items() + if not val[1].get('is_argument', False) + } + @property def value(self): - if self.env is None: - warnings.warn("requesting value from a stateless array") - return None + if self.is_stateless: + raise TypeError( + "cannot ask for the value of a stateless object " + "(use _get_value() from the top-level object instead)") if self.name not in self.env.keys(): return None return self.env[self.name] @value.setter def value(self, val): + if self.is_stateless: + raise TypeError("cannot set the value of a stateless object") self.env[self.name] = val def _to_repr_dict(self): return {'class': self.__class__, 'name': self.name, 'expr': self.expr, - 'arguments': self.arguments, - 'bound_arguments': self.bound_arguments, - 'intermediaries': self.intermediaries, + 'namespace': self.namespace, 'env': self.env} # }}} End base lazy obj class @@ -245,14 +449,13 @@ class Array(LazyObjectBase): .. attribute:: shape - overall shape of array expression. When ndim is given, - the shape itself can be a tuple of :class:`Unsigned`, - otherwise an :class:`Array` with proper assumptions and - integral dtype. + overall shape of array expression. + Tuple of :class:`Unsigned` with length equal to ``ndim``. .. attribute:: inames index names for iterating over the array. + Tuple of :class:`pymbolic.Variable`'s. .. attribute:: domain_expr @@ -270,303 +473,184 @@ class Array(LazyObjectBase): a :class:`bool` indicating whether array elements are integers - .. attribute:: assumptions + .. attribute:: integer_domain - assumptions represented as an isl set, well-defined only if - is_integral is True + an ISL set if is_integral is True, otherwise None. The set represents + the domain for the full array. It is used to express conditions like + non-negative integers, multiples of two, etc. """ _counter = 0 _name_prefix = '__lappy_array_' + # {{{ constructor + def __init__(self, name=None, shape=None, domain_expr=None, **kwargs): - # default expr + """Constructor. + """ if 'expr' not in kwargs: - kwargs['expr'] = None + kwargs['expr'] = None # will make a default expr if name is not None: kwargs['name'] = name + + stateless = kwargs.pop('stateless', False) + if not stateless: + # name default env + if 'env' not in kwargs or kwargs['env'] is None: + kwargs['env'] = default_env() + super(Array, self).__init__(**kwargs) - self._ndim = to_lappy_unsigned(kwargs.pop("ndim", None)) self.domain_expr = domain_expr + ndim = kwargs.pop("ndim", None) - if self.ndim is None: + if ndim is None and shape is not None: # infer ndim from given shape - # NOTE: use fixed ndim only shape = to_lappy_shape(shape) - self._ndim = to_lappy_unsigned(int(len(shape))) - - self.inames = kwargs.pop("inames", self._make_default_inames()) - - if self.expr is None: - self.expr = self._make_default_expr() - - if shape is None: - self._shape = self._make_default_shape() - else: - self._shape = to_lappy_shape(shape) - - self.dtype = kwargs.pop("dtype", None) - self.is_integral = kwargs.pop("is_integral", False) - self.assumptions = kwargs.pop("assumptions", None) + ndim = int(len(shape)) - # scalars - if self.shape == () or self.ndim == 0: - assert self.shape == () and self.ndim == 0 + if ndim is None: + raise ValueError("ndim cannot be determined") - # arrays with concrete ndim try: - ndim_int = int(self.ndim) - assert ndim_int == len(self.shape) + ndim = int(ndim) except ValueError: - pass - - # integer arrays - if self.assumptions is not None: - assert self.is_integral - - def _to_repr_dict(self): - repr_dict = super(Array, self)._to_repr_dict() - repr_dict.update({ - 'ndim': self.ndim, - 'inames': self.inames, - 'domain': self.domain_expr, - 'shape': self.shape, - 'dtype': self.dtype, - 'is_integral': self.is_integral, - 'assumptions': self.assumptions}) - return repr_dict + raise ValueError('ndim must be a known non-negative integer') - @property - def ndim(self): - if self._ndim == 0: - return 0 + if ndim > 0: + self.inames = kwargs.pop("inames", self._make_default_inames(ndim)) else: - return self._ndim.value - - @property - def shape(self): - """Returns the shape as a tuple of ints if the value is known, or strings of - names otherwise. - """ - return tuple(s.value if s.value is not None else s.name for s in self._shape) + self.inames = () + assert isinstance(self.inames, tuple) - @property - def size(self): - return self.shadow_size() + if self.expr is None: + self.expr = self._make_default_expr(ndim) - def shadow_size(self, axes=None): - """When axes is None, computes the full size. Otherwise computes - the size of the projected array on the given axes. - """ - if axes is None: - axes = list(range(self.ndim)) + if shape is None: + lapshape = self._make_default_shape(ndim) + for s in lapshape: + self._decl(s) + self.shape = tuple(s.name for s in lapshape) else: - # so that generators can be understood - axes = list(axes) - shape_vars = tuple( - var(self._shape[ax].name) if isinstance(self._shape[ax], Array) - else str(self._shape[ax]) - for ax in axes) - size_expr = Product(shape_vars) - shape_ctx = {s.name: s.value for s in self._shape - if isinstance(s, Array) and s.value is not None} - try: - return evaluate(size_expr, shape_ctx) - except UnknownVariableError: - return size_expr - - @property - def is_leaf(self): - """Return True if the array expression is a flat-out read. - """ - return self.expr == self._make_default_expr() + if all(isinstance(s, str) for s in shape): + self.shape = tuple(shape) + for s in self.shape: + # make new symbols if not given by the namespace arg + if s not in self.namespace: + sap = to_lappy_unsigned(s) + self._decl(sap) + else: + lapshape = to_lappy_shape(shape) + self.shape = tuple(s.name for s in lapshape) + for s in lapshape: + self._meld_namespace(s.namespace, s.as_stateless()) + self._meld_env(s.env) - def reshape(self, newshape, order='C', name=None, inames=None): - """Reshape. - """ - from lappy.core.basic_ops import reshape - return reshape( - self, newshape=newshape, - order=order, name=name, inames=inames) + self.dtype = kwargs.pop("dtype", None) + self.is_integral = kwargs.pop("is_integral", False) + if self.is_integral: + self.integer_domain = kwargs.pop("integer_domain", None) + else: + self.integer_domain = None - def transpose(self, axes=None, name=None): - """Transpose. - """ - from lappy.core.basic_ops import transpose - return transpose(self, axes, name) + # }}} End constructor - def __iter__(self): - """Iterator support if ndim and shape are all fixed. - Its behavior mimics numpy.ndarray.__iter__() - """ - if isinstance(self.ndim, Array): - raise ValueError("Cannot iterate an array with vaiable ndim: %s)" - % str(self.ndim)) - if isinstance(self.shape, Array): - raise ValueError("Cannot iterate an array with vaiable shape: %s" - % str(self.shape)) - if isinstance(self.shape[0], Array): - raise ValueError("Cannot iterate an array with vaiable shape: %s " - "(where the first component is a variable)" - % str(self.shape)) - - sym_indices = self._make_axis_indices() - for idx in range(self.shape[0]): - if self.assumptions is None: - cp_assumptions = None - else: - cp_assumptions = self.assumptions.copy() - # FIXME: project assumptions & arguments - yield self.__class__( - name=self.name + str(list([idx])), - ndim=self.ndim - 1, - shape=self.shape[1:], dtype=self.dtype, - expr=Substitution(self.expr, sym_indices[0], idx), - is_integral=self.is_integral, - arguments=self.arguments.copy(), - bound_arguments=self.bound_arguments.copy(), - env=self.env.copy(), - preconditions=list(self.preconditions), - assumptions=cp_assumptions, - value=self.value[idx], - ) + # {{{ copy constructors, with_xxx(), as_xxx() - def nditer(self): - """Flat iterator support if ndim and shape are all fixed. - Its behavior mimics numpy.ndarray.nditer() - """ - if isinstance(self.ndim, Array): - raise ValueError("Cannot iterate an array with vaiable ndim: %s)" - % str(self.ndim)) - if isinstance(self.shape, Array): - raise ValueError("Cannot iterate an array with vaiable shape: %s" - % str(self.shape)) - if any([isinstance(lenaxis, Array) for lenaxis in self.shape]): - raise ValueError("Cannot iterate an array with vaiable shape: %s " - "(contains variable components)" - % str(self.shape)) - - from pytools import indices_in_shape - sym_indices = self._make_axis_indices() - for ind in indices_in_shape(self.shape): - if self.assumptions is None: - cp_assumptions = None - else: - cp_assumptions = self.assumptions.copy() - # FIXME: project assumptions & arguments - yield self.__class__( - name=self.name + str(ind), - ndim=0, shape=(), dtype=self.dtype, - expr=Substitution(self.expr, sym_indices, ind), - is_integral=self.is_integral, - arguments=self.arguments.copy(), - bound_arguments=self.bound_arguments.copy(), - env=self.env.copy(), - preconditions=list(self.preconditions), - assumptions=cp_assumptions, - value=self.value[ind], - ) + def __copy__(self): + """Used if copy.copy is called on an array. Returns a copy of the array. - def _make_default_inames(self): - """Make the default inames, simply as lookup expressions. + Equivalent to ``a.copy(stateless=False)``. """ - return self._make_axis_indices() + return self.copy(stateless=False) - def _make_default_shape(self): - """Make a default shape based of name and ndim + def __deepcopy__(self): + """Used if copy.deepcopy is called on an array. """ - # zero-dimensional, aka a scalar - if self.ndim == 0: - return () - - # symbolic-dimensions - if isinstance(self.ndim, Unsigned): - assert self.ndim.ndim == 0 - return Lookup(Variable(self.name), "shape") - - # constant positive dimensions - assert self.ndim > 0 - names = ['__' + self.name + '_shape_%d' % d for d in range(self.ndim)] - return tuple( - Unsigned( - name=names[d], - expr=var(names[d])) - for d in range(self.ndim)) + raise NotImplementedError() def copy(self, stateless=False): """Returns a new copy of self, where the attributes undergo a shallow copy (expect for the value). + + Note: the copy operation does not make copies of the captured array data. """ - if self.assumptions is None: - cp_assumptions = None + if self.integer_domain is None: + cp_integer_domain = None else: - cp_assumptions = self.assumptions.copy() + cp_integer_domain = self.integer_domain.copy() - if stateless: + if stateless or self.is_stateless: + stateless = True env = None else: env = self.env.copy() + nms = self.namespace.copy() + return self.__class__( - name=self.name, ndim=self._ndim, + name=self.name, ndim=self.ndim, inames=self.inames, - shape=self._shape, dtype=self.dtype, + shape=self.shape, + dtype=self.dtype, + stateless=stateless, expr=self.expr, is_integral=self.is_integral, domain_expr=self.domain_expr, - arguments=self.arguments.copy(), - bound_arguments=self.bound_arguments.copy(), - intermediaries=self.intermediaries.copy(), env=env, + namespace=nms, preconditions=list(self.preconditions), - assumptions=cp_assumptions, + integer_domain=cp_integer_domain, ) + def with_state(self, namespace=None, env=None): + """Swap out the captured state. + """ + sarr = self.as_stateless() + sarr.env = env + sarr.namespace = namespace + return sarr + def as_stateless(self): """Returns a stateless copy of self. """ return self.copy(stateless=True) - def with_name(self, name, rename_arguments=False): + def with_name(self, name, rename_arguments=True): """Rename the array object. Returns a new (shallow) copy of self with the new name. + + :param rename_arguments: If true, all occurrences of the old name will + be replaced. Otherwise, the new array takes the old array as an + arguments and its value is assigned to the new name. """ new_arr = self.copy() if name == self.name or name is None: - warnings.warn("name is unchanged, just making a copy") + # warnings.warn("name is unchanged, just making a copy") return new_arr - if name in self.arguments: - raise ValueError("Name %s is already taken by an argument" % name) - elif name in self.bound_arguments: - raise ValueError("Name %s is already taken by a (bound) argument" % name) - elif name in self.intermediaries: - raise ValueError("Name %s is already taken by an intermediary" % name) + if name in self.namespace: + raise ValueError("Name %s is already taken" % name) if name in self.env: raise ValueError( "Name %s is already taken by a captured variable" % name) - if self.name in self.env: + if (self.name in self.env) and rename_arguments: new_arr.env[name] = new_arr.env.pop(self.name) - # when the argument is self-referencing + # new self-reference + new_arr.namespace[name] = new_arr.namespace.pop(self.name) + + # old self-reference if rename_arguments: - if self.name in self.arguments: - new_arr.arguments[self.name] = \ - new_arr.arguments.pop(self.name) - elif self.name in self.bound_arguments: - new_arr.bound_arguments[self.name] = \ - new_arr.bound_arguments.pop(self.name) new_arr.expr = pymbolic.substitute(new_arr.expr, {var(self.name): var(name)}) else: - if self.name in self.arguments: - new_arr.arguments.pop(self.name) - new_arr.arguments[self.name] = self.as_stateless() - elif self.name in self.bound_arguments: - new_arr.bound_arguments.pop(self.name) - new_arr.bound_arguments[self.name] = self.as_stateless() + new_arr.namespace[self.name] = ( + self.as_stateless(), self.namespace[self.name][1]) + new_arr.namespace[name] = ( + None, self.namespace[self.name][1].copy()) + new_arr.namespace[name][1]['is_argument'] = False new_arr.name = name return new_arr @@ -579,8 +663,11 @@ class Array(LazyObjectBase): assert len(inames) == len(self.inames) for name in inames: - if self.has_name(name): - raise ValueError("name %s is already taken") + if self._has_name(name.name) and (name not in self.inames): + raise ValueError( + "name %s is already taken" + % name + ) iname_map = { iname: new_iname @@ -590,539 +677,851 @@ class Array(LazyObjectBase): new_arr.expr = pymbolic.substitute(self.expr, iname_map) return new_arr - def has_name(self, name): - """Check whether a given name is taken by one of the following: - - - self.name - - self.inames - - self.arguments - - self.bound_arguments - - self.intermediaries - - self.env - """ - return False - - def _shape_names(self): - """Shape information, stringified to be legal names. - """ - for s in self._shape: - assert isinstance(s, Unsigned) - return tuple(s.name for s in self._shape) - - @property - def _shape_str(self): - """A shape description string used for loopy kernel construction. - """ - shape_names = self._shape_names() - return ', '.join(shape_names) - - @property - def _shape_expr(self): - """Values for the shape (if not None), otherwise expressions. - """ - return tuple( - s.value if s.value is not None else s.expr - for s in self._shape) - - def check_preconditions(self, context): - """Call checkers in the list of preconditions. - - :param context: a dict which is passed to all checkers. - """ - failed_checks = [] - for checker in self.preconditions: - try: - res = checker(context) - if res is not None: - if isinstance(res, bool): - if res: - pass - else: - raise PreconditionNotMetError( - "precondition checker failed: %s" % str(checker)) - else: - raise ValueError( - "cannot understand the return value " - "of the precondition checker %s" % str(checker)) - except Exception as e: # noqa: W0703 - print(e) - failed_checks.append(checker) - - err_msgs = [] - for fc in failed_checks: - msg = ("Precondition %s not met, which was imposed at\n\n" - % str(checker)) - if hasattr(fc, "frame"): - msg = msg + '\n'.join(format_list(extract_stack( - fc.frame))) - err_msgs.append(msg) - if len(failed_checks) > 0: - precon_err_msg = ( - "%d out of %d preconditions are not met" - % (len(failed_checks), len(self.preconditions))) - for im, msg in enumerate(err_msgs): - precon_err_msg = (precon_err_msg + '\n\n' - + '[%d / %d]: ' % (im + 1, len(failed_checks)) + msg) - raise PreconditionNotMetError(precon_err_msg) - def with_dtype(self, dtype, new_name=None): """Returns a copy of self with the specified dtype. """ raise NotImplementedError() - def with_shape_data(self, shape_data, new_name=None): - """Returns a copy of the array with concrete shape. Assuming concrete ndim. - ndim is unchanged by assigning the shape data, meaning that the inputs, when - not None, must be positive. - - :param shape_data: a tuple of shape data. Use None to skip setting certain - axes. - - NOTE: if the shape has nontrivial expression (not just a reference to a - name), lappy will not solve the equation for the "independent variables". - The behavior in such cases is undefined. + def with_shape_data(self, shape_data, + new_name=None, allow_complex_expr=False): + """Returns a copy of the array with concrete shape. ``ndim`` is + unchanged by assigning the shape data, meaning that the inputs, + when not None, must be positive. + + :param shape_data: a tuple of shape data. Use None to skip setting + certain axes. + + NOTE: if the shape has nontrivial expression (not just a reference + to a name), lappy will not solve the equation for the "independent + variables". The behavior in such cases is undefined. For example, if + the shape has an expression ``(n * m) % p``, this method will assign + a value for the whole expression, and the contraint among ``m, n, p`` + given by this expression is ignored and may cause unpredictable + issues. """ - if not isinstance(shape_data, (list, tuple)): - raise ValueError() - if self.ndim != len(shape_data): - raise ValueError() + assert isinstance(shape_data, (list, tuple)) + assert self.ndim == len(shape_data) + + def is_trivial(expr): + # TODO: better ways of spotting nontrivial exprs? + return expr.__class__.init_arg_names == ('name',) - for s in self._shape: - if False: # TODO: test for nontrivial exprs - warnings.warn( - "the shape variable %s of %s has nontrivial expression, " - " and setting its value may yield undefined behavior" - % (s.name, self.name)) + assert allow_complex_expr or all( + is_trivial(self._extract(s).expr) for s in self.shape) if new_name is None: new_name = 'bound%d_' % self.__class__._counter + self.name - new_arr = self.with_name(new_name) + new_arr = self.with_name(new_name, False) - delta_dict = {} - for s, s_val in zip(new_arr._shape, shape_data): - if s_val is None: - continue + for s, sv in zip(self.shape, shape_data): + if s in new_arr.env and new_arr.env[s] != sv: + raise ValueError('found pre-existing shape data') else: - assert isinstance(s_val, int) and s_val > 0 - delta_dict[s.name] = (s.value, s_val) - s.value = s_val - - # iterate through arguments and intermediaries and - # propagate the assignments to all symbols of the same name - # - # (they may not refer to the same object, but may have the same name) - # - for arg in new_arr.arguments.values(): - if arg is None: - continue - for sym_s in arg._shape: - if sym_s.name in delta_dict: - if sym_s.value == delta_dict[sym_s.name][0]: # old value - sym_s.value = delta_dict[sym_s.name][1] - else: - # the value is updated, since it points to the - # same Unsigned obj - assert sym_s.value == delta_dict[sym_s.name][1] - - for imd in new_arr.intermediaries.values(): - for sym_s in imd._shape: - if sym_s.name in delta_dict: - if sym_s.value == delta_dict[sym_s.name][0]: # old value - sym_s.value = delta_dict[sym_s.name][1] - else: - # the value is updated, since it points to the - # same Unsigned obj - assert sym_s.value == delta_dict[sym_s.name][1] + new_arr.env[s] = sv return new_arr - def with_data(self, new_name=None, **kwargs): + def with_data(self, **kwargs): """Binds specific data to a slots in the arguments, resulting in a new (shallow) copy of self with the bound arguments catupred into its environment. """ - if new_name is None: - new_name = 'bound%d_' % self.__class__._counter + self.name + if 'new_name' in kwargs: + raise ValueError("cannot rename while binding data") for key in kwargs.keys(): - if key not in self.arguments.keys(): - if key in self.bound_arguments.keys(): - # bounding cannot be undone - # (since un-inference is not possible) - raise ValueError("argument named %s is already bound" % key) + if key not in self.namespace: + raise ValueError("cannot bind an undefined name: %s" % key) - raise ValueError("argument named %s is not accepted" % key) + if key in self.env.keys(): + # bounding cannot be undone + # (since un-inference is not possible) + raise ValueError("argument named %s is already bound" + % key) + if not self.namespace[key][1].get('is_argument', False): + raise ValueError("trying to bind data to a non-argument %s" + % key) + + # make a copy new_arr = self.copy() - new_arr.name = new_name for key, val in kwargs.items(): - new_arr.bound_arguments[key] = new_arr.arguments.pop(key) + if hasattr(val, 'dtype'): if key == self.name: new_arr.dtype = val.dtype else: - new_arr.bound_arguments[key].dtype = val.dtype + new_arr.namespace[key][0].dtype = val.dtype + if hasattr(val, 'shape'): - if key == self.name: - new_arr._set_shape_values(val.shape) - else: - new_arr.bound_arguments[key]._set_shape_values(val.shape) + new_arr._set_shape_values(val.shape, key) + if hasattr(val, 'ndim'): if key == self.name: assert new_arr.ndim == val.ndim - assert len(new_arr._shape) == val.ndim + assert len(new_arr.shape) == val.ndim else: - assert new_arr.bound_arguments[key].ndim == val.ndim + assert new_arr.namespace[key][0].ndim == val.ndim assert len(val.shape) == val.ndim - new_arr.env[key] = val - - # rename value of self - new_arr.env[new_name] = new_arr.env.pop(self.name, None) + new_arr.env[key] = val # actual data return new_arr - def get_data_mapping(self, knl=None): - """Make a data mapping using known data, tailored for giving inputs to - the loopy kernel. Returns all known information if knl is None. - - :param knl: the loopy kernel - """ - data_map = {} - shapeval_expansion_list = [] - - # gather captured data - for arr_name, varr in self.env.items(): - if arr_name in ['isl_ctx', 'cl_ctx', 'cu_ctx']: - continue - - if arr_name == self.name: - # only specify output shape, and let the backend to do the malloc - for out_s in self._shape: - if out_s.value is None: - shapeval_expansion_list.append(out_s) - data_map[out_s.name] = out_s.value - - if isinstance(varr, Array) and varr.value is not None: - data_map[arr_name] = varr.value - elif isinstance(varr, np.ndarray): - data_map[arr_name] = varr - elif varr is None: - pass # self - else: - raise RuntimeError("unrecogonized captured variable %s" % arr_name) - - # try to get as much extra data that loopy wants as possible - # also evaluates the shape expressions - if knl is None: - # gather all known shapes - for s in self._shape: - if s.name in data_map: - continue - if s.value is not None: - data_map[s.name] = s.value - for arg in chain(self.arguments.values(), - self.bound_arguments.values(), self.intermediaries.values()): - if arg is None: - continue # skip self - assert isinstance(arg, Array) - for s in arg._shape: - if s.value is None: - shapeval_expansion_list.append(s) - data_map[s.name] = s.value - else: - for argname, arg in knl.arg_dict.items(): - if argname in data_map: - continue - if argname in self._shape_names(): - shape_val = self._shape[self._shape_names().index(argname)] - if shape_val.value is None: - shapeval_expansion_list.append(shape_val) - data_map[argname] = shape_val.value - for arr_arg in self.bound_arguments.values(): - if arr_arg is None: - continue - if argname in arr_arg._shape_names(): - shape_val = arr_arg._shape[ - arr_arg._shape_names().index(argname)] - if shape_val.value is None: - shapeval_expansion_list.append(shape_val) - data_map[argname] = shape_val.value - for arr_imd in self.intermediaries.values(): - if argname in arr_imd._shape_names(): - shape_val = arr_imd._shape[ - arr_imd._shape_names().index(argname)] - if shape_val.value is None: - shapeval_expansion_list.append(shape_val) - data_map[argname] = shape_val.value - # no need to search arguments - assert len(self.arguments) == 0 - - # evaluate shape expressions - for se in shapeval_expansion_list: - try: - seval = evaluate(se.expr, data_map) - except UnknownVariableError: - warnings.warn( - "cannot get value for %s prior to calling the loopy kernel" - % se.name) - seval = None - data_map[se.name] = seval - - # purge None-valued entries - entries_to_purge = [] - for key, val in data_map.items(): - if val is None: - entries_to_purge.append(key) - for key in entries_to_purge: - data_map.pop(key) - return data_map - - def eval(self, shape_dtype=np.int32): - """Evaluates the array expression after binding all the arguments. - Note that this function has side effects: it sets self.value to the - evaluation result. If the results is available, it is returned without - re-computing. - """ - if self.value is not None: - logging.info("cache hit when evaluating %s" % self.name) - return self.value + # }}} End copy constructors, with_xxx(), as_xxx() - if len(self.arguments) > 0: - raise ValueError("cannot evaluate a partially-bound expression " - "(missing data for %s)" % ', '.join(self.arguments.keys())) - - # Step 1: - # identify codegen statements (loopy assignments) - # and raise exception if ther are multiple statements - # FIXME - n_insns = 1 - if n_insns > 1: - raise TypeError("the array %s cannot be evaluated without temporaries " - "(use the compiler instead) " % self.name) - - # Step 2: - # make the initial loop domain by multiplying the index space of - # the array with hidden axes (e.g. in the case of inner products) - pass - - # Step 3: - # make the conditional expression that accounts for - # - # - sliced LHS A[1:n-1:2, :] = ... - # - # - conditional assignment A[A>0] = ... - # - # Note that this is done on the RHS, e.g. - # A[i, j] = if((i >= 1) and (i < n - 1) and ((i - 1) % 2 == 0), ..., ...) - # A[i, j] = if(A[i, j] > 0, ..., ...) - # FIXME - rhs_expr = self.expr - - # Step 4 - # Account for fancy indexing on the LHS A[B] = ... - # by adding a separate loop domain over the index space of B - # and making corresponding LHS expressions - # FIXME - lhs_expr = var(self.name)[self.inames] - loop_domain = self._generate_index_domains() - - # Step 5 - # Make and run loopy kernel - # - # FIXME: collect all assumptions and pass to loopy - # - kernel_args = [] - - # array args - for arr_name, varr in self.env.items(): - if arr_name in ['isl_ctx', 'cl_ctx', 'cu_ctx']: - continue - if arr_name == self.name: - continue # handle self later - elif arr_name in self.arguments.keys(): - if self.arguments[arr_name] is None: - pass # self - else: - assert isinstance(self.arguments[arr_name], Array) - kernel_args.append(lp.GlobalArg(arr_name, - shape=self.arguments[arr_name]._shape_str)) - elif arr_name in self.bound_arguments.keys(): - if self.bound_arguments[arr_name] is None: - pass # self - else: - assert isinstance(self.bound_arguments[arr_name], Array) - kernel_args.append(lp.GlobalArg(arr_name, - shape=self.bound_arguments[arr_name]._shape_str)) - else: - raise RuntimeError("cannot find shape information of %s" % arr_name) - - # output array args - kernel_args.append(lp.GlobalArg(self.name, shape=self._shape_str)) - - # let loopy do the inference - kernel_args.append('...') - - knl = lp.make_kernel( - loop_domain, - [ - lp.Assignment(lhs_expr, rhs_expr) - ], - kernel_args, - lang_version=(2018, 2)) - - knl = lp.set_options(knl, return_dict=True) - - # help with dtype inference - extra_dtype_dict = {} - for argname, arg in knl.arg_dict.items(): - if arg.dtype is None: - # array args - if argname == self.name: - if self.dtype is not None: - extra_dtype_dict[argname] = self.dtype - continue - - if argname in self.arguments: - if self.arguments[argname].dtype is not None: - extra_dtype_dict[argname] = self.arguments[argname].dtype - continue - elif argname in self.bound_arguments: - if self.bound_arguments[argname].dtype is not None: - extra_dtype_dict[argname] = \ - self.bound_arguments[argname].dtype - continue - - # shape args - if argname in self._shape_names(): - extra_dtype_dict[argname] = shape_dtype - for arr_arg in self.bound_arguments.values(): - if arr_arg is None: - continue - if argname in arr_arg._shape_names(): - extra_dtype_dict[argname] = shape_dtype - - knl = lp.add_and_infer_dtypes(knl, extra_dtype_dict) - data_map = self.get_data_mapping(knl) - - # TODO: to test everything that the kernel sees, - # let loopy infer more context first (e.g. shape of input arrays) - self.check_preconditions(context=data_map) - - queue = cl.CommandQueue(self.env['cl_ctx']) - evt, lp_res = knl(queue, **data_map) - - self.value = lp_res[self.name] - return self.value - - def _make_default_expr(self): - """The default expression is the array it self. - """ - sym_arr = var(self.name) - sym_indices = self._make_axis_indices() - if self.ndim > 0: - expr = Subscript(sym_arr, index=sym_indices) - else: - # scalar - expr = sym_arr - return expr + # {{{ public (numpy-compliant) api - def _set_shape_values(self, shape): - """Instantiate shape with concrete values. + @property + def ndim(self): + return len(self.shape) - This method should not be called from a user. Users should use - Array.with_shape_data() instead. - """ - for s in tuple(shape): - assert int(s) == s - assert s >= 0 - for s, sym_s in zip(shape, self._shape): - assert isinstance(sym_s, Unsigned) - sym_s.value = s + @property + def size(self): + return self._shadow_size() - def __add__(self, other): - if not isinstance(other, Array): - warnings.warn( - "Implicit conversion of %s to Lappy array" - % str(type(other))) - other_arr = to_lappy_array(other) + def itemsize(self): + """TODO: dtype expression + """ + if self.dtype is None: + # use numpy defaults + return np.dtype(self.dtype).itemsize else: - other_arr = other - from lappy.core.ufuncs import add - return add(self, other_arr) + return np.dtype(self.dtype).itemsize - __radd__ = __add__ - - @staticmethod - def _generic_constructor(obj, data=None): - if data is None: - data = obj - shape = list(obj.shape) - name = "__array_" + str(Array._counter) - _, indices = Array._create_indices(shape) - expr = var(name)[tuple(indices)] - data_mapping = {name: data} - shapes = {name: (shape, 0)} - global_index = Array._init_indices(shape, "index") - original_index = {name: indices} - Array._counter += 1 - new_array = {"expr": expr, "data_mapping": data_mapping, - "shape": shape, "shapes": shapes, - "global_index": global_index, - "original_index": original_index} - return new_array - - @staticmethod - def _generate_oned_index(new_shape): - '''Generates a 1D index for the array based off the new shape''' - temp_idx = [] - _, reshape_indices = Array._create_indices(new_shape) - for i in range(len(new_shape) - 1): - idx = reshape_indices[i] * functools.reduce(mul, new_shape[i + 1:], 1) - temp_idx.append(idx) - temp_idx.append(reshape_indices[len(reshape_indices) - 1]) - oned_idx = sum(temp_idx) - return oned_idx - - @staticmethod - def _generate_reshape_index(oned_idx, child_index, self_shape): - '''Generates the indices for the reshaped array using the 1D index''' - diff_shape = len(self_shape) - len(child_index) - if diff_shape < 0: - diff_shape = 0 - new_indices = [] - oned_index = pymbolic.substitute(oned_idx, {}) - for i in range(diff_shape): - block_size = functools.reduce(mul, self_shape[i + 1:], 1) - index = oned_index // block_size - oned_index = oned_index - (index * block_size) - for i in range(len(child_index) - 1): - block_size = functools.reduce(mul, self_shape[i + diff_shape + 1:], 1) - index = oned_index // block_size - if child_index[i] == 0: - new_indices.append(0) - else: - new_indices.append(index) - oned_index = oned_index - (index * block_size) - if child_index[len(child_index) - 1] == 0: - new_indices.append(0) - else: - new_indices.append(oned_index % self_shape[len(self_shape) - 1]) - return new_indices + def nbytes(self): + return self.size * self.itemsize - def _generate_index_domains(self): - """Make the loop domain for looping over the array's index space. - For this method to work, ndim must be concrete. - """ - assert np.isscalar(self.ndim) and abs(self.ndim - int(self.ndim)) == 0 + def __len__(self, order='C'): + if self.ndim == 0: + return 0 - inames = tuple(ina.name for ina in self.inames) - shape_names = self._shape_names() + if order == 'C': + return self.shape[0] + elif order == 'F': + return self.shape[-1] + else: + raise ValueError( + "order must be either 'C' or 'F' (got %s)" % str(order)) - assert len(inames) == self.ndim - assert len(shape_names) == self.ndim + # {{{ basic ops + + def reshape(self, newshape, order='C', name=None, inames=None): + """Reshape. + """ + from lappy.core.basic_ops import reshape + return reshape( + self, newshape=newshape, + order=order, name=name, inames=inames) + + def transpose(self, axes=None, name=None): + """Transpose. + """ + from lappy.core.basic_ops import transpose + return transpose(self, axes, name) + + T = transpose + + def swapaxes(self, axis1, axis2): + """Return a view of the array with axis1 and axis2 interchanged. + """ + raise NotImplementedError() + + def squeeze(self): + """Remove single-dimensional entries from the shape of a. + """ + raise NotImplementedError() + + def repeat(self): + """Repeat elements of an array. + """ + raise NotImplementedError() + + def fill(self): + raise NotImplementedError() + + def astype(self): + raise NotImplementedError() + + def resize(self, new_shape, order='C'): + """Change shape and size of array. + Shrinking: array is flattened, resized, and reshaped. + Enlarging: as above, with missing entries filled with zeros. + """ + raise NotImplementedError() + + def flatten(self): + raise NotImplementedError() + + ravel = flatten # unlike number, the two are equivalent in lappy + + # }}} End basic ops + + # {{{ data io + + def dump(self, with_env=True): + """Dump the full lazy array object into a file. + """ + raise NotImplementedError() + + def dumps(self, with_env=True): + """Dump the full lazy array object into a string. + """ + raise NotImplementedError() + + # }}} End data io + + # {{{ getter/setter + + def item(self, *args): + """Copy an element of an array to a standard Python scalar + and return it. + """ + raise NotImplementedError() + + def itemset(self, *args): + """Insert scalar into an array + (scalar is cast to array’s dtype, if possible). + """ + raise NotImplementedError() + + def take(self): + """Return an array formed from the elements of a at the + given indices. + """ + raise NotImplementedError() + + def put(self, indices, values, mode='raise'): + """Set a.flat[n] = values[n] for all n in indices. + """ + raise NotImplementedError() + + def __getitem__(self, indices): + """Returns a new array by sub-indexing the current array. + """ + if isinstance(indices, str): + # structured dtype support + raise NotImplementedError() + + if not isinstance(indices, tuple): + # advanced indexing with an array + assert isinstance(indices, Array) + assert indices.ndim == self.ndim + indices = (indices, ) + else: + # indexing / mixture of indexing and masking + assert isinstance(indices, tuple) + + from lappy.core.basic_ops import SubArrayFactory + arr_fac = SubArrayFactory(self, indices) + sub_arr = arr_fac() + return sub_arr + + def __setitem__(self, index, value): + """Set item described by index. + """ + raise NotImplementedError() + + def __delitem__(self, key): + """Delete self[key]. + """ + raise NotImplementedError() + + def __iter__(self): + """Iterator support if ndim and shape are all fixed. + Its behavior mimics numpy.ndarray.__iter__() + """ + raise NotImplementedError() + + def flat(self): + """Flat iterator support if ndim and shape are all fixed. + Its behavior mimics numpy.ndarray.nditer() + """ + raise NotImplementedError() + + # }}} End getter/setter + + # {{{ numpy protocols + + def __array__(self): + """(For numpy’s dispatch). When converting to a numpy array with + ``numpy.array`` or ``numpy.asarray``, the lazy array is evaluated + and the result is returned. + """ + res = self.eval() + if isinstance(res, np.ndarray): + return res + else: + return res.get() + + @property + def __array_interface__(self): + """The array interface (array protocol) for data buffer access. + """ + if self.is_stateless: + raise ValueError("array value is unknown") + else: + raise NotImplementedError() + + @classmethod + def __array_function__(cls, func, types, args, kwargs): + """Array function protocol. + https://numpy.org/neps/nep-0018-array-function-protocol.html + """ + raise NotImplementedError() + + @classmethod + def __array_ufunc__(cls, ufunc, method, *inputs, **kwargs): + """Array ufunc protocol. + https://docs.scipy.org/doc/numpy/reference/ufuncs.html + """ + raise NotImplementedError() + + # }}} End numpy protocols + + # {{{ arithmetics + + def __add__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import add + return add(self, other_arr) + + def __radd__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import add + return add(other_arr, self) + + def __sub__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import subtract + return subtract(self, other_arr) + + def __rsub__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import subtract + return subtract(other_arr, self) + + def __mul__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import multiply + return multiply(self, other_arr) + + def __rmul__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import multiply + return multiply(other_arr, self) + + def __div__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import divide + return divide(self, other_arr) + + __truediv__ = __div__ + + def __rtruediv__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import divide + return divide(other_arr, self) + + def __divmod__(self, value): + raise NotImplementedError() + + def __rdivmod__(self, value): + raise NotImplementedError() + + def __pow__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import power + return power(self, other_arr) + + def __rpow__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import power + return power(other_arr, self) + + def __floordiv__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import floordiv + return floordiv(self, other_arr) + + def __rfloordiv__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import floordiv + return floordiv(other_arr, self) + + def __mod__(self, other): + raise NotImplementedError() + + def __rmod__(self, other): + raise NotImplementedError() + + def __iadd__(self, value): + raise NotImplementedError() + + def __isub__(self, value): + raise NotImplementedError() + + def __imul__(self, value): + raise NotImplementedError() + + def __idiv__(self, value): + raise NotImplementedError() + + def __itruediv__(self, value): + raise NotImplementedError() + + def __ipow__(self, value): + raise NotImplementedError() + + def __ifloordiv__(self, value): + raise NotImplementedError() + + def __imod__(self, value): + raise NotImplementedError() + + def __matmul__(self, value): + raise NotImplementedError() + + def __rmatmul__(self, value): + raise NotImplementedError() + + def __imatmul__(self, value): + raise NotImplementedError() + + def __lshift__(self): + raise NotImplementedError() + + def __rlshift__(self): + raise NotImplementedError() + + def __ilshift__(self): + raise NotImplementedError() + + def __rshift__(self): + raise NotImplementedError() + + def __rrshift__(self): + raise NotImplementedError() + + def __irshift__(self): + raise NotImplementedError() + + # }}} End arithmetics + + # {{{ comparisons + + def __eq__(self, other): + """Test for equality. Returns ``True`` if two Arrays have + the same name. Returns a new Array with comparison expression + otherwise. + """ + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + + if self.name == other_arr.name: + return True + else: + from lappy.core.ufuncs import equal + return equal(self, other_arr) + + def __ne__(self, other): + return not self == other + + def __le__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import less_equal + return less_equal(self, other_arr) + + def __lt__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import less + return less(self, other_arr) + + def __ge__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import greater_equal + return greater_equal(self, other_arr) + + def __gt__(self, other): + if not isinstance(other, Array): + other_arr = to_lappy_array(other) + else: + other_arr = other + from lappy.core.ufuncs import greater + return greater(self, other_arr) + + # }}} End comparisons + + # {{{ logical ops + + def __and__(self, value): + raise NotImplementedError() + + def __rand__(self, value): + raise NotImplementedError() + + def __or__(self, value): + raise NotImplementedError() + + def __ror__(self, value): + raise NotImplementedError() + + def __xor__(self, value): + raise NotImplementedError() + + def __rxor__(self, value): + raise NotImplementedError() + + def __iand__(self, value): + raise NotImplementedError() + + def __ior__(self, value): + raise NotImplementedError() + + def __ixor__(self, value): + raise NotImplementedError() + + def __contains__(self, key): + raise NotImplementedError() + + def __invert__(self): + raise NotImplementedError() + + def all(self, axis=None, out=None, keepdims=False): + raise NotImplementedError() + + def any(self, axis=None, out=None, keepdims=False): + raise NotImplementedError() + + # }}} End logical ops + + # {{{ other maths + + def __abs__(self): + """|self| + """ + from lappy.core.ufuncs import absolute + return absolute(self) + + def __pos__(self): + """+self + """ + raise NotImplementedError() + + def __neg__(self): + """-self + """ + from lappy.core.ufuncs import negative + return negative(self) + + def __index__(self): + """Called to implement operator.index(), and whenever Python needs to + losslessly convert the numeric object to an integer object (such as in + slicing, or in the built-in bin(), hex() and oct() functions). + """ + if self.is_integral: + raise NotImplementedError() + else: + raise ValueError() + + def __complex__(self): + """Convert size-1 arrays to Python scalars. + """ + raise NotImplementedError() + + def __float__(self): + """Convert size-1 arrays to Python scalars. + """ + raise NotImplementedError() + + def __int__(self): + """Convert size-1 arrays to Python scalars. + """ + raise NotImplementedError() + + def round(self): + """Return a with each element rounded to the given number of + decimals. + """ + raise NotImplementedError() + + def conjugate(self): + raise NotImplementedError() + + conj = conjugate + + def real(self): + raise NotImplementedError() + + def imag(self): + raise NotImplementedError() + + def max(self, axis=None, out=None, keepdims=False, initial=None, where=True): + raise NotImplementedError() + + def min(self): + raise NotImplementedError() + + def sort(self): + raise NotImplementedError() + + def searchsorted(self): + """Find indices where elements of v should be inserted in a + to maintain order. + """ + raise NotImplementedError() + + def mean(self): + raise NotImplementedError() + + def sum(self): + raise NotImplementedError() + + def var(self): + raise NotImplementedError() + + def std(self): + raise NotImplementedError() + + def nonzero(self): + raise NotImplementedError() + + def partition(self): + raise NotImplementedError() + + def prod(self): + """Return the product of the array elements over the given axis. + """ + raise NotImplementedError() + + product = prod + + def ptp(self): + """Peak to peak (maximum - minimum) value along a given axis. + """ + raise NotImplementedError() + + def anom(self): + """Compute the anomalies (deviations from the arithmetic mean) + along the given axis. + """ + raise NotImplementedError() + + def argmax(self): + raise NotImplementedError() + + def argmin(self): + raise NotImplementedError() + + def argpartition(self): + raise NotImplementedError() + + def argsort(self): + raise NotImplementedError() + + def clip(self): + """Clip (limit) the values in an array. + """ + raise NotImplementedError() + + def cumprod(self): + """Return the cumulative product of elements along a given axis. + """ + raise NotImplementedError() + + def cumsum(self): + """Return the cumulative sum of the elements along the given axis. + """ + raise NotImplementedError() + + def diagonal(self): + raise NotImplementedError() + + def trace(self): + raise NotImplementedError() + + def dot(self): + raise NotImplementedError() + + # }}} End other maths + + # }}} End public (numpy-compliant) api + + # {{{ extended public api + + def pin(self, name=None, cse=False, insert_barrier=False): + """Declare that the array shall be treated as a temporary. + """ + if name is None: + name = self.name + + if cse: + self.namespace[name][1]['is_cse'] = True + + if insert_barrier: + self.namespace[name][1]['has_barrier'] = True + + self.namespace[name][1]['is_temporary'] = True + + # }}} End extended public api + + # {{{ private api + + def _has_trivial_expr(self): + """Whether the expression is more complex than a leaf node. + """ + # An array + if isinstance(self.expr, Subscript): + if self.expr.aggregate.__class__.init_arg_names == ('name',): + for idx in self.expr.index_tuple: + if self.expr.aggregate.__class__.init_arg_names == ('name',): + pass + else: + return False + return True + + # A scalar + return self.expr.__class__.init_arg_names == ('name',) + + def _get_shape_vals(self, axis=None): + """Returns the shape values where available, None otherwise. + """ + return tuple(self.env[s.name] if s.name in self.env else None + for s in self.shape) + + def _set_shape_values(self, shape, name=None): + """Instantiate shape of a namespace member with concrete values. + """ + for s in tuple(shape): + assert int(s) == s + assert s >= 0 + + if name is None or name == self.name: + sym_shape = self.shape + elif name in self.namespace: + sym_shape = self._extract(name).shape + else: + raise ValueError() + + for s, sym_s_name in zip(shape, sym_shape): + sym_s = self._extract(sym_s_name) + if sym_s.name in self.env: + if self.env[sym_s_name] != s: + warnings.warn( + "captured shape value for %s is overwritten " + "(%s --> %s)" % ( + sym_s.name, + str(self.env[sym_s.name]), str(s))) + self.env[sym_s_name] = s + + def _get_value(self, name=None): + """Try to get the value of a captured variable. Returns None + if the process fails. + """ + if name is None: + return self.value + + if self.is_stateless: + return None + + if name in self.env: + return self.env[name] + else: + return None + + def _get_shape_vals_or_strs(self): + """Returns the shape as a tuple of ints if the value is known, or + strings of names otherwise. + """ + return tuple( + self.env[s] if s in self.env else s for s in self.shape) + + def _make_default_expr(self, *args): + """The default expression is the array it self. + + :param ndim: optional, used when the array shape is not set. + """ + if hasattr(self, 'shape'): + ndim = self.ndim + else: + ndim, = args + + sym_arr = var(self.name) + if ndim > 0: + # array + sym_indices = self._make_axis_indices(*args) + expr = Subscript(sym_arr, index=sym_indices) + else: + # scalar + expr = sym_arr + return expr + + def _generate_index_domains(self): + """Make the loop domain for looping over the array's index space. + """ + if self.ndim == 0: + # make a set with a single element for scalar output + return [ + isl.BasicSet.read_from_str( + isl.DEFAULT_CONTEXT, + "{ [%s] : %s >= 0 and %s < 1 }" + % ((self.name + '_dummy_iname', ) * 3)) + ] + + inames = tuple(ina.name for ina in self.inames) + shape_names = self.shape + + assert len(inames) == self.ndim + assert len(shape_names) == self.ndim v = isl.make_zero_and_vars(inames, shape_names) @@ -1133,190 +1532,116 @@ class Array(LazyObjectBase): & v[inames[i]].lt_set(v[shape_names[i]])) return s.get_basic_sets() - @staticmethod - def _transpose_list_args(indices, args): - '''Transposes the axes of the array based off the given args''' - new_indices = [] - for i in range(len(indices)): - if args[i] >= len(indices): - raise AssertionError("Given argument does not match with shape axes") - new_indices.append(indices[args[i]]) - return new_indices - - @staticmethod - def _add_shapes(obj1, obj2): - '''Adds two shapes of arrays together, following the rules of broadcasting''' - shape = list(obj1) - min_dims = min(len(obj1), len(obj2)) - first_idx = len(obj1) - 1 - second_idx = len(obj2) - 1 - for j in range(min_dims): - # If the dimensions are equal or one of the dimensions is 1, then - # broadcasting is legal. - first_shape_dim = shape[first_idx] - second_shape_dim = obj2[second_idx] - - if ( - first_shape_dim == 1 - or second_shape_dim == 1 - or (first_shape_dim == second_shape_dim)): - if first_idx <= (len(shape) - 1): - shape[first_idx] = max(first_shape_dim, second_shape_dim) - else: - shape.insert(0, max(first_shape_dim, second_shape_dim)) - else: - return None - first_idx -= 1 - second_idx -= 1 - # take care of the rest of the dims once we have filled min_dims - if len(shape) < len(obj2): - for k in range(second_idx, -1, -1): - shape.insert(0, obj2[k]) - return shape - - @staticmethod - def _add_dicts(obj1, obj2): - '''Combines two dictionaries''' - data_map = {} - for k, v in obj1.items(): - data_map[k] = v - for k, v in obj2.items(): - data_map[k] = v - return data_map - - def _make_axis_indices(self): + def _make_axis_indices(self, *args): """Make the inames as expressions based of array name and ndim + + :param ndim: optional, used when the array shape is not set yet. """ # zero-dimensional, aka a scalar - if self.ndim == 0: + if hasattr(self, 'shape'): + ndim = self.ndim + else: + ndim, = args + + if ndim == 0: return () # ndim is None (up to further inference) - if self.ndim is None: - # return Lookup(Variable(self.name), "inames") - return var('__%s_inames' % self.name) - - # symbolic-dimensions - if isinstance(self.ndim, Unsigned): - assert self.ndim.ndim == 0 + if ndim is None: # return Lookup(Variable(self.name), "inames") return var('__%s_inames' % self.name) # constant positive dimensions - assert self.ndim > 0 + assert ndim > 0 # shape = Lookup(Variable(self.name), "inames") # return tuple(Subscript(shape, d) for d in range(self.ndim)) return tuple(var('__%s_inames_%d' % (self.name, d)) - for d in range(self.ndim)) - - @staticmethod - def _create_indices(shape): - ''' - Generates standard indices for an array based off the shape. - Note that a 0 is used for axes of shape 1. - Indexing starts from the "back" of the array. - ''' - names = [] - indices = [] - idx = 0 - for i in range(len(shape) - 1, -1, -1): - name = "index_" + str(idx) - if shape[i] == 1: - indices.insert(0, 0) - else: - indices.insert(0, var(name)) - names.insert(0, name) - idx += 1 - return names, indices - - @staticmethod - def _init_indices(set, type): - ''' - Initializes a set of indices for the array, based on the given set and type: - "index" for array expression or "s" for shape indices - i.e. [index_2, index_1, index_0] or [s2, s1, s0] - Note that this function ignores using 0 for axes of shape 1. - ''' - indices = [] - idx = 0 - for i in range(len(set) - 1, -1, -1): - name = type + "_" + str(idx) - indices.insert(0, var(name)) - idx += 1 - return indices - - @staticmethod - def _check_data(data_mapping): - '''Returns True if the array expression has data for all arrays''' - for k, v in data_mapping.items(): - if v is None: - return False - return True + for d in range(ndim)) - @staticmethod - def _check_shape(shape): - '''Returns True if the shape has defined values for all indices''' - for i in range(len(shape)): - if not isinstance(shape[i], int): - return False - return True + def _to_repr_dict(self): + repr_dict = super(Array, self)._to_repr_dict() + repr_dict.update({ + 'ndim': self.ndim, + 'inames': self.inames, + 'domain': self.domain_expr, + 'shape': self.shape, + 'dtype': self.dtype, + 'is_integral': self.is_integral, + 'integer_domain': self.integer_domain}) + return repr_dict - @staticmethod - def _update_shape(shapes): - ''' - Updates the shape based off the shapes of each array in the - array expression - ''' - new_shape = [] - for k, v in shapes.items(): - if v is None: - return None - if v[1] == 0: - shape = v[0] - else: - shape = v[1] - new_shape = Array._add_shapes(new_shape, shape) - if new_shape is None: - return None - return new_shape - - @staticmethod - def _transpose_indices(global_index, args): - ''' - Creates a dictionary mapping between the global index and the corresponding - global index in the transposed array - ''' - transpose_dict = {} - for i in range(len(args)): - transpose_dict[global_index[i]] = global_index[args.index(i)] - return transpose_dict - - @staticmethod - def _generate_args(global_index): - ''' - Generates the args representing the permutation of the axes for the - global index. Permutation of axes is defined as any variation from the - standard indices originally assigned - ''' - args = [] - for i in range(len(global_index) - 1, -1, -1): - name = "index_" + str(i) - args.append(global_index.index(var(name))) - return args - - @staticmethod - def _update_zero_indices(old_index, new_shape): - ''' - Generates a dictionary mapping between the old index of the array - and the new shape, using 0's for shapes of 1 - ''' - expr_dict = {} - for i in range(len(old_index)): - if new_shape[i] == 1: - expr_dict[old_index[i]] = 0 - else: - expr_dict[old_index[i]] = old_index[i] - return expr_dict + def _shadow_size(self, axes=None): + """When axes is None, computes the full size. Otherwise computes + the size of the projected array on the given axes. + + If the shasow size of not computable, returns an expression for it. + """ + if axes is None: + axes = list(range(self.ndim)) + else: + # so that generators can be understood + axes = list(axes) + shape_vars = tuple(var(self.shape[ax]) for ax in axes) + size_expr = Product(shape_vars) + try: + return evaluate(size_expr, self.env) + except UnknownVariableError: + return size_expr + + @property + def _is_leaf(self): + """Return True if the array expression is a flat-out read. + """ + return self.expr == self._make_default_expr() + + def _make_default_inames(self, *args): + """Make the default inames, simply as lookup expressions. + """ + return self._make_axis_indices(*args) + + def _make_default_shape(self, ndim): + """Make a default shape based of name and ndim. + Returns the shape tuple of ``Unsigned``. + """ + # zero-dimensional, aka a scalar + if ndim == 0: + return () + + # symbolic-dimensions + if isinstance(ndim, Unsigned): + assert ndim.ndim == 0 + return Lookup(Variable(self.name), "shape") + + # constant positive dimensions + assert ndim > 0 + if len(self.name) > 1 and self.name[:2] == '__': + name_prefix = '' + else: + name_prefix = '__' + names = [name_prefix + self.name + '_shape_%d' % d + for d in range(ndim)] + return tuple( + Unsigned( + name=names[d], + expr=var(names[d]), + env=None, namespace=None) + for d in range(ndim)) + + @property + def _shape_str(self): + """A shape description string used for loopy kernel construction. + """ + return ', '.join(self.shape) + + @property + def _shape_expr(self): + """Values for the shape (if not None), otherwise expressions. + """ + return tuple( + self.env[s] if s in self.env else self._extract(s).expr + for s in self.shape) + + # }}} End private api # }}} End array class @@ -1325,18 +1650,32 @@ class Array(LazyObjectBase): # for more readable names and handy constructors -class Int(Array): - """Integer (index type). +class Scalar(Array): + """Scalar. """ _counter = 0 - _name_prefix = '__lappy_int_' + _name_prefix = '__lappy_scalar_' def __init__(self, name=None, **kwargs): if name is not None: kwargs['name'] = name - kwargs['shape'] = () kwargs['ndim'] = 0 + super(Scalar, self).__init__(**kwargs) + + def __getitem__(self, indices): + raise TypeError("cannot sub-index a scalar") + + +class Int(Scalar): + """Integer (index type). + """ + _counter = 0 + _name_prefix = '__lappy_int_' + + def __init__(self, name=None, **kwargs): + if name is not None: + kwargs['name'] = name super(Int, self).__init__(**kwargs) self.is_integral = True @@ -1355,8 +1694,8 @@ class Unsigned(Int): if name is not None: kwargs['name'] = name super(Unsigned, self).__init__(**kwargs) - self.assumptions = isl.BasicSet.read_from_str( - self.env['isl_ctx'], + self.integer_domain = isl.BasicSet.read_from_str( + isl.DEFAULT_CONTEXT, "{ [%s] : %s >= 0 }" % ((self.name, ) * 2)) assert self.is_integral @@ -1380,7 +1719,7 @@ def to_lappy_array(arr_like, name=None, base_env=None): if np.isscalar(arr_like): try: - dtype = np.dtype(arr_like) + dtype = np.dtype(type(arr_like)).type except TypeError: dtype = None @@ -1388,11 +1727,11 @@ def to_lappy_array(arr_like, name=None, base_env=None): if int(arr_like) == arr_like: arr_class = Int else: - arr_class = Array + arr_class = Scalar - arr = arr_class(name=name, dtype=dtype, env=base_env) + arr = arr_class(name=name, dtype=dtype, env=base_env, value=arr_like) arr.env[arr.name] = arr_like - arr.arguments = dict() + return arr if isinstance(arr_like, np.ndarray): @@ -1416,44 +1755,105 @@ def to_lappy_array(arr_like, name=None, base_env=None): raise ValueError("Cannot convert the input to a Lappy Array.") -def to_lappy_unsigned(unsigned_like, name=None, base_env=None): +def to_lappy_scalar(scalar_like, name=None, base_env=None, dtype=np.float64): + """Do nothing if the array is already an :class:`Scalar`. + Make an :class:`Scalar` with captured values for an actual number. + """ + if base_env is None: + base_env = default_env() + + if scalar_like is None: + return Scalar(name=name, env=base_env, dtype=dtype) + + if np.isscalar(scalar_like): + try: + scalar_like = float(scalar_like) + return Scalar(value=scalar_like, env=base_env, dtype=dtype) + except ValueError: + # a string of name, not value + assert isinstance(scalar_like, str) + return Scalar(name=scalar_like, env=base_env, dtype=dtype) + + if isinstance(scalar_like, Array): + assert scalar_like.ndim == 0 + assert scalar_like.shape == () + + if not isinstance(scalar_like, Scalar): + scalar_like = Scalar( + name=scalar_like.name, value=scalar_like.value, + dtype=dtype, + env=scalar_like.env.copy(), + namespace=scalar_like.namespace.copy()) + + return scalar_like + + if isinstance(scalar_like, pymbolic.primitives.Expression): + return Scalar(expr=scalar_like, dtype=dtype, env=base_env) + + raise ValueError("Cannot convert the input to a Lappy Scalar.") + + +def to_lappy_unsigned(unsigned_like, name=None, base_env=None, dtype=np.int32): """Do nothing if the array is already an :class:`Unsigned`. Make an :class:`Unsigned` with captured values for an actual number. - - Specifically, returns 0 if the value is 0. This is to avoid infinite - recursion when using Array type to represent shape and ndim. """ if base_env is None: base_env = default_env() + if name is None: + name = make_default_name(Unsigned) + if unsigned_like is None: - return Unsigned(name=name, env=base_env) + return Unsigned( + name=name, env=base_env, + namespace={name: (None, {'is_argument': True})}, + dtype=dtype) if np.isscalar(unsigned_like): - if unsigned_like == 0: - return 0 + try: + unsigned_like = int(unsigned_like) + except ValueError: + # a string of name, not value + return Unsigned( + name=unsigned_like, env=base_env, + namespace={ + unsigned_like: (None, {'is_argument': True})}, + dtype=dtype) + if abs(int(unsigned_like)) - unsigned_like == 0: - lunsigned = Unsigned(name=name, value=unsigned_like, env=base_env) + lunsigned = Unsigned( + name=name, value=unsigned_like, env=base_env, + namespace={name: (None, {'is_argument': True})}, + dtype=dtype) lunsigned.env[lunsigned.name] = unsigned_like return lunsigned else: - raise ValueError("%s is not a non-negative integer" % str(unsigned_like)) + raise ValueError( + "%s is not a non-negative integer" % str(unsigned_like)) if isinstance(unsigned_like, Array): assert unsigned_like.ndim == 0 assert unsigned_like.shape == () assert unsigned_like.is_integral - assert unsigned_like.value is not None if not isinstance(unsigned_like, Unsigned): unsigned_like = Unsigned( - name=unsigned_like.name, value=unsigned_like.value, - env=unsigned_like.env.copy()) + name=unsigned_like.name, + value=unsigned_like.value, + dtype=dtype, + env=unsigned_like.env.copy(), + namespace=unsigned_like.namespace.copy(), + expr=unsigned_like.expr, + arguments=unsigned_like.arguments.copy(), + bound_arguments=unsigned_like.bound_arguments.copy(), + temporaries=unsigned_like.temporaries.copy(), + preconditions=unsigned_like.preconditions + ) return unsigned_like if isinstance(unsigned_like, pymbolic.primitives.Expression): - return Unsigned(expr=unsigned_like) + return Unsigned(expr=unsigned_like, dtype=dtype, env=base_env) raise ValueError("Cannot convert the input to a Lappy Unsigned.") @@ -1470,8 +1870,6 @@ def to_lappy_shape(shape): Unsigned objects. - Comma/space-separated string. - Nonnegativity assumption is added to the input where the shape is - parametrized. """ if shape is None: return tuple() @@ -1488,9 +1886,9 @@ def to_lappy_shape(shape): int_c = int(c) if not int_c == abs(int_c): raise ValueError("cannot use %s as a shape parameter" % str(c)) - components[i] = Unsigned(name=None, value=int(c)) + components[i] = to_lappy_unsigned(int(c)) except ValueError: - components[i] = Unsigned(c) + components[i] = to_lappy_unsigned(c) return tuple(components) assert isinstance(shape, (tuple, list)) @@ -1498,25 +1896,41 @@ def to_lappy_shape(shape): for i, c in enumerate(shape): if isinstance(c, Array): components[i] = c - if not isinstance(c, Unsigned): - warnings.warn( - "Implicit conversion of %s to unsigned (non-negative) " - "integer" - % str(c)) else: - try: - int_c = int(c) - if not int_c == abs(int_c): - raise ValueError("cannot use %s as a shape parameter" % str(c)) - components[i] = Unsigned(name=None, value=int(c)) - except ValueError: - assert isinstance(c, str) - components[i] = Unsigned(c) - except TypeError: - assert isinstance(c, Expression) - components[i] = Unsigned(expr=c) + components[i] = to_lappy_unsigned(c, name=None) return tuple(components) +def isscalar(obj): + """Like np.isscalar, but also works for lazy objects. + """ + if np.isscalar(obj): + return True + else: + return isinstance(obj, Int) + # }}} End digest shape specifiers + +# {{{ misc utils + + +def get_internal_name_prefix(name): + """If the name is already an internal name (starting with '__'), return ''; + otherwise return '__'. + """ + if len(name) > 1 and name[:2] == '__': + name_prefix = '' + else: + name_prefix = '__' + return name_prefix + + +def make_default_name(obj_class): + """Make a default name for the lazy object class. + """ + name = '%s%d' % (obj_class._name_prefix, obj_class._counter) + obj_class._counter += 1 + return name + +# }}} End misc utils diff --git a/lappy/core/basic_ops.py b/lappy/core/basic_ops.py index 038546a86bad446606526bda249bf5bb0199fef3..4e4083bf89fa4e09663e29c699469e5eae9296a1 100644 --- a/lappy/core/basic_ops.py +++ b/lappy/core/basic_ops.py @@ -21,14 +21,354 @@ THE SOFTWARE. """ import numpy as np from pymbolic import var, evaluate, substitute -from pymbolic.primitives import Variable, Product -from numpy import AxisError +from pymbolic.primitives import Variable, Product, Subscript -from lappy.core.array import to_lappy_shape, to_lappy_unsigned +from lappy.exceptions import AxisError +from lappy.core.ufuncs import minimum, maximum +from lappy.core.conditional import conditional +from lappy.core.array import ( + Array, to_lappy_shape, to_lappy_unsigned, isscalar) from lappy.core.preconditions import ( EvaluationPrecondition, make_size_conservation_condition) +# {{{ masking + + +class MaskedArrayFactory(object): + """Factory class of the sub-arrays produced by indexing with a boolean + array with equal shape. + + Filtered arrays cannot be created as stand-alone expressions; they + are used solely for the lhs of an assignment. For example, to execute + + A[A > 0] = B[A > 0] + + It is achieved by overloading __setitem__() of the LHS. + + If the user asks for the value of A[A > 0], the kernel will compute + the full array A and the boolean array (A > 0). And then lappy will + simply use numpy to give the results. + """ + + def __init__(self, base, mask): + assert isinstance(base, Array) + assert isinstance(mask, Array) + assert base.ndim == mask.ndim + self.base_arr = base + self.mask_arr = mask + # TODO: precondition on equal shapes + + def __call__(self): + # TODO: make a special kind of array, whose evaluation returns multiple + # results, which form the final reults after some post-processing. + # + # e.g., A[A > 0] + raise NotImplementedError() + +# }}} End masking + +# {{{ indexing + + +class SubArrayFactory(object): + """Factory class of the sub-arrays produced by indexing with a tuple of + (generalized) indices. + + For the syntax and semantics are the same as numpy arrays + https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html + + .. attribute:: base_arr + + the base array being indexed. + + .. attribute:: selectors + + a list of selection objects. + + .. attribute:: ndim + + ndim of the output array. + + .. attribute:: advanced_indexing + + list of bools, whether is using advanced indexing for each axis. + + .. attribute:: boolean_masking + + list of bools, whether is masked with a boolean array for each axis. + Meaningful only if advanced_indexing is True. + """ + + # {{{ constructor + + def __init__(self, base, selobj=None): + """ + :param base_arr: the base array. + :param selectors: the selection object. + """ + assert isinstance(base, Array) + self.base_arr = base + + # normalize the selectors list + if isinstance(selobj, tuple): + selist = list(selobj) + elif isinstance(selobj, (int, Array, list, np.ndarray)): + selist = [selobj, ] + elif selobj is None: + selist = [] + else: + raise TypeError() + + assert isinstance(selist, list) + n_ell = 0 + for idx in selist: + if idx is Ellipsis: + n_ell += 1 + + if n_ell > 1: + raise IndexError( + "an index can only have a single ellipsis ('...')") + elif n_ell == 0: + # hand un-addressed axes with a trailing ellipsis + selist.append(Ellipsis) + n_ell += 1 + + self.selectors = selist + self.ndim = self._get_out_ndim() + self._expand_ellipsis() + + self.advanced_indexing = self._prepare_advanced_indexing() + assert len(self.advanced_indexing) == self.ndim + + self.boolean_masking = self._prepare_boolean_masking() + assert len(self.boolean_masking) == self.ndim + + # }}} End constructor + + # {{{ consume the selection object + + def _get_out_ndim(self): + """Compute the rank (ndim) of the output array. + """ + ndim = 0 + n_annihil_dim = 0 + for sel in self.selectors: + if isinstance(sel, slice): # slicing + ndim += 1 + elif sel is None: # np.newaxis + ndim += 1 + elif isscalar(sel): # scalar indexing + n_annihil_dim += 1 + elif isinstance(sel, (tuple, list, np.ndarray, Array)): + sel_arr = np.array(sel) + if sel_arr.dtype is np.dtype(bool) or sel_arr.dtype is bool: + # boolean indexing "flattens" + ndim += 1 + else: + # integer indexing keeps the dimension + ndim += sel_arr.ndim + else: + if sel is not Ellipsis: + raise TypeError() + + return max(ndim, self.base_arr.ndim, n_annihil_dim) + + def _prepare_advanced_indexing(self): + """Identify if the selection object is/contains advanced indexing + for each output axis. + """ + advanced = [False, ] * self.ndim + for isel, sel in enumerate(self.selectors): + if isinstance(sel, Array): + advanced[isel] = True + return advanced + + def _prepare_boolean_masking(self): + """Identify if the selection object is/contains advanced indexing + via boolean masking for each output axis. + """ + bool_masked = [False, ] * self.ndim + for isel, sel in enumerate(self.selectors): + if isinstance(sel, Array): + if sel.dtype == bool or self.dtype == np.bool_: + bool_masked[isel] = True + return bool_masked + + def _complete_slice(self, slc, iaxis): + """Fill in all information for a partially specified slice. + + :param slc: an (incomplete) slice object. + :param iaxis: the index of the axis being sliced. + """ + assert isinstance(slc, slice) + dim_name = self.base_arr.shape[iaxis] + dim = self.base_arr._extract(dim_name) + + step_none = slc.step is None + if step_none: + step = 1 + else: + step = slc.step + + start = slc.start + start_none = slc.start is None + if not start_none: + start = conditional(start < 0, start + dim, start) + + stop = slc.stop + stop_none = slc.stop is None + if not stop_none: + stop = conditional(stop < 0, stop + dim, stop) + + start_p = 0 if start_none else maximum(0, minimum(dim, start)) + stop_p = dim if stop_none else maximum(start_p, minimum(dim, stop)) + + start_m = dim - 1 if start_none \ + else maximum(-1, minimum(dim - 1, start)) + stop_m = -1 if stop_none else maximum(-1, minimum(start_m, stop)) + + start = conditional(step > 0, start_p, start_m) + stop = conditional(step > 0, stop_p, stop_m) + + return slice(start, stop, step) + + def _expand_ellipsis(self): + """Expand ellipsis into proper number of slices. + """ + ell_id = self.selectors.index(Ellipsis) + fullslice = slice(None, None, None) + self.selectors[ell_id:ell_id + 1] = ( + (fullslice, ) * (self.base_arr.ndim - len(self.selectors) + 1)) + + # }}} End consume the selection object + + # {{{ basic indexing + + def _basic_getitem(self, basic_selectors): + """Tthe basic_selectors must only consist of: + slices, integers, Ellipsis, and np.newaxis (None). + """ + ca = 0 # current axis + i_newaxis = 0 + new_shape = [] + new_inames = [] + rhs_index_exprs = [] + + new_namespace = self.base_arr.namespace.copy() + new_namespace[self.base_arr.name] = ( + self.base_arr.as_stateless(), + self.base_arr.namespace[self.base_arr.name][1].copy()) + + new_name = self.base_arr._make_default_name() + + for i_sel, sel in enumerate(basic_selectors): + if sel is None: + # newaxis does not advance the ca pointer + new_shape.append(to_lappy_unsigned(1)) + new_inames.append(var('__indexing_%s_newaxis_inames_%d' + % (self.base_arr.name, i_newaxis))) + i_newaxis += 1 + + elif ca >= self.base_arr.ndim: + raise IndexError('too many indices for array') + + elif isinstance(sel, slice): + se = self._complete_slice(sel, iaxis=ca) + + se_size = (se.stop - se.start + se.step - 1) // se.step + new_shape.append(se_size) + + iname = self.base_arr.inames[ca] + new_inames.append(iname) + + if isinstance(se.start, Array): + start = se.start.expr + new_namespace[se.start.name] = ( + se.start, + se.start.namespace[se.start.name][1].copy()) + else: + start = se.start + + if isinstance(se.step, Array): + step = se.step.expr + new_namespace[se.step.name] = ( + se.step, + se.step.namespace[se.step.name][1].copy()) + else: + step = se.step + + rhs_index_exprs.append( + start + step * iname) + ca += 1 + + elif isscalar(sel): + new_shape.append(1) + iname = self.base_arr.inames[ca] + new_inames.append(iname) + rhs_index_exprs.append(sel) + ca += 1 + + else: + raise TypeError('invalid index type: %s' + % type(basic_selectors[i_sel])) + + # NOTE: since complex expressions cannot be subscrpted, + # if the expression is non trivial, this will automatically + # add a temporary. + new_expr = Subscript( + var(self.base_arr.name), + index=tuple(rhs_index_exprs)) + + obj = { + 'name': new_name, + 'inames': tuple(new_inames), + 'expr': new_expr, + 'value': None, + 'domain_expr': self.base_arr.domain_expr, + 'env': self.base_arr.env.copy(), + 'namespace': new_namespace, + 'preconditions': self.base_arr.preconditions, + 'ndim': self.ndim, + 'shape': new_shape, + 'dtype': self.base_arr.dtype, + 'is_integral': self.base_arr.is_integral, + 'integer_domain': self.base_arr.integer_domain, + } + arr = self.base_arr.__class__(**obj) + + if not self.base_arr._has_trivial_expr(): + # make previous results cse to ensure correct dependence + arr.pin(self.base_arr.name, cse=True) + + return arr + + # }}} End basic indexing + + # {{{ advanced indexing + + def _advanced_getitem(self): + """Handle advanced indexing of two cases: + + 1) array of size_t + 2) array of bool + """ + raise NotImplementedError() + + # }}} End advanced indexing + + def __call__(self): + """Create the sub-array. + """ + if any(self.advanced_indexing): + return self._advanced_getitem() + + return self._basic_getitem(self.selectors) + +# }}} End indexing + +# {{{ reshape + def reshape(array, newshape, order='C', name=None, inames=None): """Reshape an array to the given shape. @@ -95,16 +435,16 @@ def reshape(array, newshape, order='C', name=None, inames=None): """The shape must be an integer along the newaxis. """ old_size = evaluate(array.size, context) - shadow_size = evaluate(Product(shadow_expr), context) + _shadow_size = evaluate(Product(shadow_expr), context) act_new_shape = tuple( var(s.name) if hasattr(s, 'name') else s for s in newshape) - if not old_size % shadow_size == 0: + if not old_size % _shadow_size == 0: raise ValueError("cannot reshape array of size %s into shape %s" % (str(old_size), str(act_new_shape))) - return old_size % shadow_size == 0 + return old_size % _shadow_size == 0 # add pre-eval runtime checks on size consistency new_precond.append( @@ -155,40 +495,49 @@ def reshape(array, newshape, order='C', name=None, inames=None): if order == 'C': iname_subs = { - array.inames[0]: flat_id // array.shadow_size(range(1, array.ndim)) + array.inames[0]: flat_id // array._shadow_size(range(1, array.ndim)) } for iaxis in range(1, array.ndim): - flat_id = flat_id % array.shadow_size(range(iaxis, array.ndim)) - iname_subs[array.inames[iaxis]] = flat_id // array.shadow_size( + flat_id = flat_id % array._shadow_size(range(iaxis, array.ndim)) + iname_subs[array.inames[iaxis]] = flat_id // array._shadow_size( range(iaxis + 1, array.ndim)) else: assert order == 'F' iname_subs = { - array.inames[-1]: flat_id // array.shadow_size(range(array.ndim - 1)) - } + array.inames[-1]: flat_id // array._shadow_size(range(array.ndim - 1)) + } for iaxis in range(-2, -array.ndim - 1, -1): - flat_id = flat_id % array.shadow_size(range(array.ndim + iaxis + 1)) - iname_subs[array.inames[iaxis]] = flat_id // array.shadow_size( + flat_id = flat_id % array._shadow_size(range(array.ndim + iaxis + 1)) + iname_subs[array.inames[iaxis]] = flat_id // array._shadow_size( range(array.ndim + iaxis)) - new_arr = array.with_name(name) - new_arr._ndim = to_lappy_unsigned(new_ndim) + new_arr = array.with_name(name, False) new_arr.inames = inames new_arr.preconditions = new_precond - new_arr._shape = to_lappy_shape(newshape) # noqa: W0212 + newshape = to_lappy_shape(newshape) + new_arr.shape = tuple(s.name for s in newshape) + for s in newshape: + new_arr._meld_namespace(s.namespace, s.as_stateless()) + new_arr._meld_env(s.env) - if name != array.name: - # if the reference is no longer to self - if array.name in new_arr.arguments.keys(): - new_arr.arguments[array.name] = array.as_stateless() - elif array.name in new_arr.bound_arguments.keys(): - new_arr.bound_arguments[array.name] = array.as_stateless() + new_arr.namespace[array.name] = ( + array.as_stateless(), + array.namespace[array.name][1].copy()) + new_arr.namespace[name] = (None, {}) + + if not array._has_trivial_expr(): + # make previous results cse to ensure correct dependencies + new_arr.pin(array.name, cse=True) new_arr.expr = substitute(array.expr, iname_subs) return new_arr +# }}} End reshape + +# {{{ transpose + def transpose(array, axes=None, name=None): """Transpose. @@ -229,8 +578,10 @@ def transpose(array, axes=None, name=None): name = (name_prefix + array.name + '_T%d' % array.__class__._counter) # noqa: W0212 - new_arr = array.with_name(name) + new_arr = array.with_name(name, False) new_arr.inames = tuple(array.inames[axes[i]] for i in range(ndim)) - new_arr._shape = tuple(array._shape[axes[i]] for i in range(ndim)) + new_arr.shape = tuple(array.shape[axes[i]] for i in range(ndim)) return new_arr + +# }}} End transpose diff --git a/lappy/core/broadcast.py b/lappy/core/broadcast.py index 11e5fe48414b1e1e366fad1058e76ab22e805070..fda79d64f6aef9d018bc85c67414e3b61f4862d7 100644 --- a/lappy/core/broadcast.py +++ b/lappy/core/broadcast.py @@ -24,7 +24,7 @@ from functools import partial from pymbolic import evaluate, substitute, var from pymbolic.primitives import Product, Expression, Max -from lappy.core.array import Unsigned, to_lappy_shape +from lappy.core.array import to_lappy_shape, to_lappy_array, make_default_name from lappy.core.tools import is_nonnegative_int from lappy.core.preconditions import EvaluationPrecondition @@ -33,6 +33,10 @@ class BroadcastResult(object): """Encapsulates the broadcasting result. The object acts like an array when performing shape-related queries. + .. attribute:: name + + name of the broadcast, used to name the resulting arrays. + .. attribute:: base_arrays list of (lazy) arrays, whose api should support .shape and .ndim @@ -58,48 +62,64 @@ class BroadcastResult(object): size (number of elements) of the broadcast result + .. attribute:: inames + + all broadcast results share the same set of inames + .. attribute:: preconditions list of extra preconditions to make the broadcast feasible NOTE: only fixed (statically known) ndims are supported """ + _counter = 0 + _name_prefix = '__lappy_broadcast_res_' + def __init__(self, array_list): - for arr in array_list: - if not np.isscalar(arr.ndim): - raise ValueError( - "cannot broadcast %s with variable ndim" - % str(arr)) - self.base_arrays = list(array_list) + self.base_arrays = array_list self.preconditions = [] + + self.name = self._make_broadcast_name() self._broadcase_base_arrays() + self.inames = self._make_broadcast_inames() + self._make_broadcast_arrays() + def _make_broadcast_name(self): + return make_default_name(self.__class__) + + def _make_broadcast_inames(self): + if self.ndim == 0: + return () + assert self.ndim > 0 + return tuple(var('__%s_inames_%d' % (self.name, d)) + for d in range(self.ndim)) + def _broadcase_base_arrays(self): """Compute the ndim, shape and size of the broadcast result. """ ndims = [arr.ndim for arr in self.base_arrays] - self._ndim = max(ndims) + self.ndim = max(ndims) # values when available, names otherwise - shapes = [arr.shape for arr in self.base_arrays] + shapes = [arr._get_shape_vals_or_strs() for arr in self.base_arrays] # gather shape expressions expr_map = {} for arr in self.base_arrays: - for s in arr._shape: - if s.name not in expr_map: - expr_map[var(s.name)] = s.expr + for s in arr.shape: + if s not in expr_map: + expr_map[var(s)] = arr._extract(s).expr else: # same name must have the same runtime values # (may have different expressions) - def check_name_valule_consistency(context, s): + def check_name_value_consistency(context, s): s_val = evaluate(s.expr, context) another_s_val = evaluate( expr_map[var(s.name)], context) return s_val == another_s_val self.preconditions.append(EvaluationPrecondition( - partial(check_name_valule_consistency, s=s))) + partial(check_name_value_consistency, s=s))) # implicitly reshape to the same ndim # (right-aligned) @@ -190,27 +210,26 @@ class BroadcastResult(object): bshape = tuple(bshape_pre[i] for i in range(self.ndim)) self._shape_exprs = bshape - self._size_expr = Product(self._shape_exprs) + # Unsigneds of broadcast shape + self.shape = to_lappy_shape(self._shape_exprs) + def _make_broadcast_arrays(self): """After knowing ndim, shape and size symbolically, construct the broadcast copies of the input arrays. """ self.broadcast_arrays = [] - # Unsigneds of broadcast shape - lappy_shape = to_lappy_shape(self._shape_exprs) - for base_arr in self.base_arrays: dim_offset = self.ndim - base_arr.ndim assert dim_offset >= 0 # ints and strs - in_shape = base_arr.shape + in_shape = base_arr._get_shape_vals_or_strs() if dim_offset == 0: - if in_shape == self.shape: + if in_shape == self.shape_val_or_str: # unchanged self.broadcast_arrays.append(base_arr) continue @@ -223,21 +242,22 @@ class BroadcastResult(object): name = (name_prefix + name + '_broadcast%d' % base_arr.__class__._counter) # noqa: W0212 - brdc_arr = base_arr.with_name(name) + brdc_arr = base_arr.with_name(name, False) assert isinstance(self.ndim, int) and self.ndim >= 0 - brdc_arr._ndim = Unsigned(value=self.ndim) - brdc_arr._shape = lappy_shape + brdc_arr.shape = tuple(s.name for s in self.shape) + for s in self.shape: + brdc_arr._decl(s) # make new inames - brdc_arr.inames = brdc_arr._make_default_inames() + brdc_arr.inames = self.inames # make a new copy with broadcast expression # # account for semi-dynamic broadcasting, i.e., passing shape value 1 # at runtime, by fetching from (out_iname % in_shape). # TODO: when it is known that such broadcasting cannot happen, this - # may be optimized out in loopy by imposing assumptions like + # may be optimized out in loopy by imposing constraints like # out_iname < in_shape expr_mapping = {} for inm, new_inm, in_shape_i in zip( @@ -246,27 +266,21 @@ class BroadcastResult(object): (None, ) * dim_offset + base_arr._shape_expr): if inm is not None: expr_mapping[inm] = new_inm % in_shape_i + brdc_arr.expr = substitute(brdc_arr.expr, expr_mapping) - # update argument list - # if there was a self reference, the metadata needs to be captured - for argname in base_arr.arguments.keys(): - if base_arr.arguments[argname] is None: - brdc_arr.arguments[argname] = base_arr - for argname in base_arr.bound_arguments.keys(): - if base_arr.bound_arguments[argname] is None: - brdc_arr.bound_arguments[argname] = base_arr + # capture the base array + brdc_arr.namespace[base_arr.name] = ( + base_arr.as_stateless(), + base_arr.namespace[base_arr.name][1].copy() + ) self.broadcast_arrays.append(brdc_arr) assert len(self.base_arrays) == len(self.broadcast_arrays) @property - def ndim(self): - return self._ndim - - @property - def shape(self): + def shape_val_or_str(self): return tuple(s if is_nonnegative_int(s) else str(s) for s in self._shape_exprs) @@ -282,4 +296,8 @@ class BroadcastResult(object): def broadcast(*arrays): """Creates an object that mimics broadcasting. """ + arrays = list(arrays) + for aid in range(len(arrays)): + if np.isscalar(arrays[aid]): + arrays[aid] = to_lappy_array(arrays[aid]) return BroadcastResult(arrays) diff --git a/lappy/core/conditional.py b/lappy/core/conditional.py new file mode 100644 index 0000000000000000000000000000000000000000..6c79e1c17b3938fe600d6d2ca6138c4ee87cac6f --- /dev/null +++ b/lappy/core/conditional.py @@ -0,0 +1,150 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" +import numpy as np + +from pymbolic import var, substitute +from pymbolic.primitives import If +from lappy.core.array import Array, to_lappy_array +from lappy.core.tools import check_and_merge_precs, check_and_merge_envs + + +def conditional(cond, con, alt, name=None): + """If-then-else conditional expression for arrays. + + :param cond: The condition. + :param con: The consequent. + :param alt: The alternative. + :param name: (Optional) name of the resulting array. + + Behaves like the ternary operator ?: in C. + """ + # if cond is simply a known boolean, simplify the result + if isinstance(cond, (int, bool)): + if cond: + return con + else: + return alt + + if not isinstance(cond, Array): + cond = to_lappy_array(cond) + + if not isinstance(con, Array): + con = to_lappy_array(con) + + if not isinstance(alt, Array): + alt = to_lappy_array(alt) + + # con and alt must have the same ndim + assert con.ndim == alt.ndim + + # cond can be a scalar or of the same shape as con/alt + if cond.ndim > 0: + assert cond.ndim == con.ndim + + obj = dict() + + if con.domain_expr == alt.domain_expr: + obj['domain_expr'] = con.domain_expr + else: + # TODO: check/add precondition that although the + # expressions differ, the resulting domains are the same + raise NotImplementedError() + + if con.dtype is None: + if alt.dtype is None: + new_dtype = None + else: + new_dtype = alt.dtype + else: + new_dtype = np.result_type(con.dtype, alt.dtype).type + + obj['dtype'] = new_dtype + arr_class = Array + + obj['ndim'] = con.ndim + + if name is None: + name = '__if_(%s)_then_(%s)_else_(%s)' % ( + cond.name, con.name, alt.name) + + obj['name'] = name + obj['inames'] = tuple( + var('__(%s)_inames_%d' % (name, d)) + for d in range(obj['ndim'])) + + iname_maps = dict() + for iaxis in range(obj['ndim']): + iname_maps[con.inames[iaxis]] = obj['inames'][iaxis] + iname_maps[alt.inames[iaxis]] = obj['inames'][iaxis] + if cond.ndim > 0: + iname_maps[cond.inames[iaxis]] = obj['inames'][iaxis] + + # replace old inames with the same (new) ones + obj['expr'] = substitute( + If(cond.expr, con.expr, alt.expr), + iname_maps) + + # con and alt must have the same shape, cond can be of the same shape, + # or be a scalar of shape (, ). + if all(s1.name == s2.name for s1, s2 in zip(con.shape, alt.shape)): + if cond.ndim > 0: + if all(s0.name == s1.name + for s0, s1 in zip(cond.shape, con.shape)): + obj['shape'] = cond.shape + else: + # TODO: check/add preconditions as needed + raise NotImplementedError() + else: + obj['shape'] = con.shape + else: + # TODO: check/add preconditions as needed + raise NotImplementedError() + + obj['env'] = check_and_merge_envs(cond, con, alt) + obj['preconditions'] = check_and_merge_precs(cond, con, alt) + obj['namespace'] = None + + obj['is_integral'] = con.is_integral & alt.is_integral + + if not obj['is_integral']: + obj['integer_domain'] = None + else: + if con.integer_domain is None and alt.integer_domain is None: + obj['integer_domain'] = None + else: + if con.integer_domain is None: + obj['integer_domain'] = alt.integer_domain + elif alt.integer_domain is None: + obj['integer_domain'] = con.integer_domain + else: + obj['integer_domain'] = alt.integer_domain.union( + con.integer_domain) + + arr = arr_class(**obj) + + arr._meld_namespace(cond.namespace, cond.as_stateless()) + arr._meld_namespace(con.namespace, con.as_stateless()) + arr._meld_namespace(alt.namespace, alt.as_stateless()) + + return arr diff --git a/lappy/core/function_base.py b/lappy/core/function_base.py new file mode 100644 index 0000000000000000000000000000000000000000..1946e8d540fd88ac25e0f7cff791136ea7102107 --- /dev/null +++ b/lappy/core/function_base.py @@ -0,0 +1,148 @@ +from __future__ import division, absolute_import, print_function + +import numpy as np +from lappy import ndarray +from lappy.core.array import ( + default_env, + LazyObjectBase, Scalar, Unsigned, + to_lappy_scalar, to_lappy_unsigned) +from lappy.core.ufuncs import power +from lappy.core.tools import (check_and_merge_envs, check_and_merge_precs) + +__copyright__ = "Copyright (C) 2020 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__all__ = ["arange", "linspace", "logspace"] + + +def arange(start, stop=None, step=1, dtype=None, + name=None, base_env=None): + """Return evenly spaced values within a given interval. + Behaves like :func:`numpy.arange`. + + ``` + numpy.arange([start, ]stop, [step, ]dtype=None) + ``` + + :param start: Start of interval. The interval includes this value. + The default start value is 0. + :param stop: End of interval. The interval does not include this value, + :param step: Spacing between values. The default step size is 1. + :param dtype: The type of the output array. + + :param name: Name of the output. + """ + if stop is None: # the pos 1 arg is treated as stop + stop = start + start = 0 + + if base_env is None: + base_env = default_env() + + start = to_lappy_scalar(start) + stop = to_lappy_scalar(stop) + step = to_lappy_scalar(step) + + num = to_lappy_unsigned((stop - start) // step) + + obj = dict() + obj['ndim'] = 1 + obj['shape'] = (num.name, ) + obj['env'] = check_and_merge_envs(start, stop, step) + obj['namespace'] = {num.name: (num, dict())} + obj['preconditions'] = check_and_merge_precs(start, stop, step) + + if name is not None: + obj['name'] = name + + if dtype is None: + dtype = np.result_type(start, stop, step).type + obj['dtype'] = dtype + + arr = ndarray(**obj) + arr.namespace[arr.name][1]['is_argument'] = False + + iname = arr.inames[0] + arr.expr = start.expr + step.expr * iname + + arr._meld_namespace(start.namespace, start.as_stateless()) + arr._meld_namespace(stop.namespace, stop.as_stateless()) + arr._meld_namespace(step.namespace, step.as_stateless()) + + return arr + + +def linspace(start, stop, num=50, endpoint=True, retstep=False, + name=None, dtype=None): + """Return evenly spaced numbers over a specified interval. + Behaves like :func:`numpy.linspace`. + + :param start: The starting (Scalar) value. + :param stop: The end (Scalar) value. + :param num: Number of samples to generate. Default is 50. + :param endpoint: Whether the last point at ``stop`` is included. + :param retstep: If True, return the spacing as well. + + :param name: Name of the output. Auto-generate the name if None. + :param dtype: Data type of the output. + """ + if isinstance(endpoint, LazyObjectBase): + raise TypeError("boolean value 'endpoint' cannot be lazy.") + if isinstance(retstep, LazyObjectBase): + raise TypeError("boolean value 'retstep' cannot be lazy.") + if not isinstance(start, Scalar): + start = to_lappy_scalar(start) + if not isinstance(stop, Scalar): + stop = to_lappy_scalar(stop) + if not isinstance(num, Unsigned): + num = to_lappy_unsigned(num) + + if endpoint: + if num == 1: + step = 0. + else: + step = (stop - start) / (num - 1) + y = arange(0, num) * step + start + # FIXME: add this line when support indexing + # y[-1] = stop + else: + step = (stop - start) / num + y = arange(0, num) * step + start + + if retstep: + return y, step + else: + return y + + +def logspace(start, stop, num=50, endpoint=True, base=10.0): + """Return numbers spaced evenly on a log scale. + Behaves like :func:`numpy.logspace`. + + :param start: The starting (Scalar) value. + :param stop: The end (Scalar) value. + :param num: Number of samples to generate. Default is 50. + :param endpoint: Whether the last point at ``stop`` is included. + :param base: The base of the log space. Default is 10.0. + """ + y = linspace(start, stop, num=num, endpoint=endpoint) + return power(base, y) diff --git a/lappy/core/masked_array.py b/lappy/core/masked_array.py new file mode 100644 index 0000000000000000000000000000000000000000..0f51b34b4beb155dceb8c7198fb23cd9557a2222 --- /dev/null +++ b/lappy/core/masked_array.py @@ -0,0 +1,64 @@ +from __future__ import division, absolute_import, print_function + +from lappy.core.array import Array + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__doc__ = """ +============= +Masked Arrays +============= +""" + + +class MaskedArrayError(Exception): + pass + + +class MaskError(MaskedArrayError): + pass + + +class MaskedArray(Array): + """An array class with possibly masked values. + The mask array is an array of booleans that shares the shape and inames + with the (base) data array, and lives in the same closure. + + For a masked array A, A.expr and A.value means the data array. To get the + value of the final results, access A.masked_value instead. + + .. attribute:: mask_expr + + expression for the mask array + + .. attribute:: mask_value + + value for the mask array + + .. attribute:: masked_value + + value of the masked array (after post-processing) + """ + _counter = 0 + _name_prefix = '__lappy_masked_array_' + pass diff --git a/lappy/core/tools.py b/lappy/core/tools.py index 8e685b41390a9e7d199ab25d67f89bb8c206a7b8..cdb4c4180c771fe8bbe1f4679d5d767e613246de 100644 --- a/lappy/core/tools.py +++ b/lappy/core/tools.py @@ -22,17 +22,26 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ import numpy as np -from lappy.core.array import LazyObjectBase +from lappy.core.array import LazyObjectBase, Scalar +__all__ = ["check_and_merge_envs", "check_and_merge_precs", + "is_nonnegative_int", "is_constexpr", 'ScalarType'] -def check_and_merge_envs(*envs): + +ScalarType = np.ScalarType + (Scalar, ) + + +def check_and_merge_envs(*arr_list): """Merge captured environments. + Enviroments (captured lists) are dictionaries with names as keys, and + data (numpy arrays, pyopencl arrays etc.) as values. + Raises exception when there are conflicts. """ new_env = dict() - for env in envs: - for key, val in env.items(): + for arr in arr_list: + for key, val in arr.env.items(): if key in new_env: if new_env[key] is val: pass @@ -49,6 +58,20 @@ def check_and_merge_envs(*envs): return new_env +def check_and_merge_precs(*arr_list): + """Merge precondition lists. + + PLs are lists of :class:`EvaluationPrecondition` objects. + + Returns the merged argument list. + """ + mprecs = list() + for arr in arr_list: + mprecs.extend(arr.preconditions) + + return mprecs + + def is_nonnegative_int(num): if not np.isscalar(num): return False diff --git a/lappy/core/ufuncs.py b/lappy/core/ufuncs.py index 1c8d3220b898d0ddb6a8c451e3bfee82b665c5c3..2abca4c42a765b6f7b6f8803736a3fd42c3d9ca1 100644 --- a/lappy/core/ufuncs.py +++ b/lappy/core/ufuncs.py @@ -21,11 +21,14 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from pymbolic import var -from lappy.core.array import Array +import numpy as np + +from pymbolic import var, substitute +from pymbolic.primitives import Min, Max +import lappy +from lappy.core.array import Array, to_lappy_array from lappy.core.primitives import FloatingPointRemainder from lappy.core.broadcast import broadcast -from lappy.core.tools import check_and_merge_envs class UFunc(object): @@ -47,6 +50,10 @@ class UFunc(object): whether the output is integral for integral inputs. + .. attribute:: is_integral + + whether the output is always integral (regardless of inputs). + """ nargs = -1 @@ -89,10 +96,12 @@ class UFunc(object): return self.__call__(*args, **kwargs) - def __init__(self, expr_func, dtype=None, preserve_integral=False): + def __init__(self, expr_func, dtype=None, + preserve_integral=False, is_integral=False): self.dtype = dtype self.f = expr_func self.preserve_integral = preserve_integral + self.is_integral = is_integral def __call__(self, *args, **kwargs): raise NotImplementedError() @@ -104,7 +113,10 @@ class UnaryOperation(UFunc): nargs = 1 def __call__(self, a, name=None): - new_assumptions = None + new_integer_domain = None + + if not isinstance(a, Array): + a = to_lappy_array(a) # non-typed ufunc: pass-through of input's dtype if self.dtype is None: @@ -112,33 +124,34 @@ class UnaryOperation(UFunc): else: new_dtype = self.dtype - new_interm = a.intermediaries.copy() - if name in new_interm: + if name in a.namespace: raise ValueError("The name %s is ambiguous" % name) - new_interm[a.name] = a.as_stateless() + + if a.namespace[a.name][1].get('is_temporary', False): + base_expr = var(a.name) + else: + base_expr = a.expr obj = { 'name': name, 'inames': a.inames, - 'expr': self.f(a.expr), 'value': None, + 'expr': self.f(base_expr), 'value': None, 'domain_expr': a.domain_expr, - 'arguments': a.arguments.copy(), - 'bound_arguments': a.bound_arguments.copy(), 'env': a.env.copy(), + 'namespace': a.namespace.copy(), 'preconditions': list(a.preconditions), - 'intermediaries': new_interm, 'ndim': a.ndim, 'shape': a.shape, 'dtype': new_dtype, - 'is_integral': a.is_integral and self.preserve_integral, - 'assumptions': new_assumptions, + 'is_integral': + (a.is_integral and self.preserve_integral) or self.is_integral, + 'integer_domain': new_integer_domain, } - arr = Array(**obj) + arr = lappy.ndarray(**obj) - # remove self reference - if a.name in arr.arguments: - arr.arguments[a.name] = a.as_stateless() - elif a.name in arr.bound_arguments: - arr.bound_arguments[a.name] = a.as_stateless() + # update self reference + arr.namespace[a.name] = ( + a.as_stateless(), a.namespace[a.name][1].copy()) + arr.namespace[arr.name] = (None, dict()) return arr @@ -244,79 +257,67 @@ class BinaryOperation(UFunc): raise NotImplementedError() def __call__(self, a, b, name=None): + + if not isinstance(a, Array): + a = to_lappy_array(a) + + if not isinstance(b, Array): + b = to_lappy_array(b) + bres = broadcast(a, b) a, b = bres.broadcast_arrays if self.dtype is None: - # FIXME: allow numpy type promotion rules - assert a.dtype == b.dtype - new_dtype = a.dtype + if a.dtype is None: + new_dtype = b.dtype + elif b.dtype is None: + new_dtype = a.dtype + else: + new_dtype = np.result_type(a.dtype, b.dtype).type else: new_dtype = self.dtype - for ka, va in a.arguments.items(): - if ka in b.arguments.keys(): - if not va.ndim == b.arguments[ka].ndim: - raise ValueError("argument %s has different ndims in %s and %s" - % (ka, a.name, b.name)) - if not va.shape == b.arguments[ka].shape: - raise ValueError("argument %s has different shapes in %s and %s" - % (ka, a.name, b.name)) - - assert ka == va.name - assert ka == b.arguments[ka].name - - new_arglist = a.arguments.copy() - new_arglist.update(b.arguments) - - new_bound_arglist = a.bound_arguments.copy() - new_bound_arglist.update(b.bound_arguments) - - new_env = check_and_merge_envs(a.env, b.env) - new_interm = check_and_merge_envs(a.intermediaries, b.intermediaries) - - if name in new_env: - raise ValueError("The name %s is ambiguous" % name) - - new_interm[a.name] = a - new_interm[b.name] = b - - # TODO: more elegant handling of inames - assert len(a.inames) == len(b.inames) - b = b.with_inames(a.inames) - # expression handling - new_expr = self.f(a.expr, b.expr) - - new_assumptions = None + if a.namespace[a.name][1].get('is_temporary', False): + a_expr = var(a.name) + else: + a_expr = substitute( + a.expr, + {ainame: briname + for ainame, briname in zip(a.inames, bres.inames)} + ) + if b.namespace[b.name][1].get('is_temporary', False): + b_expr = var(b.name) + else: + b_expr = substitute( + b.expr, + {biname: briname + for biname, briname in zip(b.inames, bres.inames)} + ) + new_expr = self.f(a_expr, b_expr) + + new_integer_domain = None obj = { 'name': name, - 'inames': a.inames, + 'inames': bres.inames, 'expr': new_expr, 'value': None, 'domain_expr': a.domain_expr, - 'arguments': new_arglist, - 'bound_arguments': new_bound_arglist, - 'intermediaries': new_interm, - 'env': new_env, + 'namespace': None, + 'env': a.env.copy(), 'preconditions': a.preconditions + b.preconditions + bres.preconditions, 'ndim': bres.ndim, 'shape': bres._shape_exprs, 'dtype': new_dtype, - 'is_integral': all([ - a.is_integral, b.is_integral, self.preserve_integral]), - 'assumptions': new_assumptions, + 'is_integral': + all([a.is_integral, b.is_integral, self.preserve_integral]) + or self.is_integral, + 'integer_domain': new_integer_domain, } - arr = Array(**obj) - - # remove self reference - if a.name in arr.arguments: - arr.arguments[a.name] = a.as_stateless() - elif a.name in arr.bound_arguments: - arr.bound_arguments[a.name] = a.as_stateless() + arr = lappy.ndarray(**obj) - if b.name in arr.arguments: - arr.arguments[b.name] = b.as_stateless() - elif b.name in arr.bound_arguments: - arr.bound_arguments[b.name] = b.as_stateless() + arr.namespace[arr.name][1]['is_argument'] = False + arr._meld_namespace(a.namespace.copy(), a.as_stateless()) + arr._meld_namespace(b.namespace.copy(), b.as_stateless()) + arr._meld_env(b.env.copy()) return arr @@ -345,9 +346,12 @@ abs = absolute = UnaryOperation(lambda x: var('abs')(x), preserve_integral=True) fabs = UnaryOperation(lambda x: var('fabs')(x)) negative = UnaryOperation(lambda x: -x, preserve_integral=True) -floor = UnaryOperation(lambda x: var('floor')(x), preserve_integral=True) -ceil = UnaryOperation(lambda x: var('ceil')(x), preserve_integral=True) -around = UnaryOperation(lambda x: var('around')(x), preserve_integral=True) +floor = UnaryOperation(lambda x: var('floor')(x), + preserve_integral=True, is_integral=True) +ceil = UnaryOperation(lambda x: var('ceil')(x), + preserve_integral=True, is_integral=True) +around = UnaryOperation(lambda x: var('around')(x), + preserve_integral=True, is_integral=True) logical_not = UnaryOperation(lambda x: x.not_()) @@ -364,6 +368,8 @@ add = BinaryOperation(expr_func=lambda x, y: x + y, preserve_integral=True) subtract = BinaryOperation(expr_func=lambda x, y: x - y, preserve_integral=True) multiply = BinaryOperation(expr_func=lambda x, y: x * y, preserve_integral=True) divide = BinaryOperation(expr_func=lambda x, y: x / y) +floordiv = BinaryOperation(expr_func=lambda x, y: x // y, is_integral=True) +power = BinaryOperation(expr_func=lambda x, y: x**y, preserve_integral=True) equal = BinaryOperation(lambda x, y: x.eq(y)) not_equal = BinaryOperation(lambda x, y: x.ne(y)) @@ -381,6 +387,7 @@ mod = remainder = BinaryOperation(lambda x, y: x % y) # floating point mod, and returns the same sign with x fmod = BinaryOperation(lambda x, y: FloatingPointRemainder(x, y)) -# TODO: add these when needed: arctan2, alltrue, sometrue, logical_xor, hypot +minimum = BinaryOperation(lambda x, y: Min([x, y])) +maximum = BinaryOperation(lambda x, y: Max([x, y])) # }}} End binary ufuncs diff --git a/lappy/eval/__init__.py b/lappy/eval/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..536f4ffd477b10307295f0861f8bea86b6c85aee --- /dev/null +++ b/lappy/eval/__init__.py @@ -0,0 +1,57 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__doc__ = """ +Code generation and evaluation. +""" + + +class Closure(object): + """Code and data container. + """ + def __init__(self, out_names, lpknl, data_map): + self.out_names = out_names + self.kernel = lpknl + self.data_map = data_map + assert self._check_context() + + def _check_context(self): + """Check if data_map has consistent context. + """ + ctx = None + for key, val in self.data_map.items(): + if hasattr(val, 'context'): + if ctx is None: + ctx = val.context + elif ctx == val.context: + continue + else: + return False + self.context = ctx + return True + + +from lappy.eval.execution import Executor +from lappy.eval.compiler import Compiler +__all__ = ['Compiler', 'Closure', 'Executor'] diff --git a/lappy/eval/compiler.py b/lappy/eval/compiler.py new file mode 100644 index 0000000000000000000000000000000000000000..fe7d2452b07b96edb54eb9d780ad6bdb57adaad2 --- /dev/null +++ b/lappy/eval/compiler.py @@ -0,0 +1,456 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2020 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import warnings +import functools +from traceback import format_list, extract_stack + +import numpy as np +import loopy as lp +import pyopencl as cl + +from pymbolic import var, substitute, evaluate +from pymbolic.mapper import CombineMapper +from pymbolic.mapper.evaluator import UnknownVariableError + +from lappy.core.array import Array +from lappy.eval import Closure + + +class PreconditionNotMetError(Exception): + pass + +# {{{ expression mappers + + +class TemporaryCollector(CombineMapper): + """Find temporaries in the expression. + Note that it is possible that a temporary appears multiple + times within different inames. + + .. attribute:: arr + + reference to the lazy array object. + """ + def __init__(self, arr): + self.arr = arr + + def combine(self, values): + import operator + return functools.reduce(operator.add, values, list()) + + def map_constant(self, expr): + if hasattr(expr, 'name'): + if expr.name in self.arr.namespace: + tags = self.arr.namespace[expr.name][1] + is_temp = tags.get('is_temporary', False) + if is_temp: + return [expr.name, ] + + return list() + + map_algebraic_leaf = map_constant + +# }}} End expression mappers + + +# {{{ compiler + + +class Compiler(object): + def __init__(self, check_preconditions=False): + self.check_preconditions = check_preconditions + self._counter = 0 + + def __call__(self, *args): + """When given a few of :class:`Array` objects, make one loopy + kernel that evaluates all of them. + """ + if len(args) > 1: + raise NotImplementedError() + + out_names = [] + for arr in args: + assert isinstance(arr, Array) + out_names.append(arr.name) + + knl, context = self.compile(*args) + data_map = self._get_data_mapping(arr, knl) + context.update(data_map) + + # FIXME: using with_name breaks precondition checking. + # Possible Solution: keep track of name changes + if self.check_preconditions: + self._check_preconditions(arr.preconditions, context) + + return Closure(out_names, knl, data_map) + + def _check_preconditions(self, preconditions, context): + """Call checkers in the list of preconditions. + + :param context: a dict which is passed to all checkers. + """ + failed_checks = [] + for checker in preconditions: + try: + res = checker(context) + if res is not None: + if isinstance(res, bool): + if res: + pass + else: + raise PreconditionNotMetError( + "precondition checker failed: %s" + % str(checker)) + else: + raise ValueError( + "cannot understand the return value " + "of the precondition checker %s" + % str(checker)) + except Exception as e: # noqa: W0703 + print(e) + failed_checks.append(checker) + + err_msgs = [] + for fc in failed_checks: + msg = ("Precondition %s not met, which was imposed at\n\n" + % str(checker)) + if hasattr(fc, "frame"): + msg = msg + '\n'.join(format_list(extract_stack( + fc.frame))) + err_msgs.append(msg) + if len(failed_checks) > 0: + precon_err_msg = ( + "%d out of %d preconditions are not met" + % (len(failed_checks), len(preconditions))) + for im, msg in enumerate(err_msgs): + precon_err_msg = (precon_err_msg + '\n\n' + + '[%d / %d]: ' % (im + 1, len(failed_checks)) + msg) + raise PreconditionNotMetError(precon_err_msg) + + def _get_data_mapping(self, arr, knl=None): + """Make a data mapping using known data, tailored for giving inputs to + the loopy kernel. Returns all known information if knl is None. + + :param knl: the loopy kernel + """ + data_map = {} + expansion_list = [] + + # gather captured data + for arr_name, varr in arr.env.items(): + + if arr_name == arr.name: + # only specify output shape, and let the backend to do the malloc + for out_s in arr.shape: + if out_s not in arr.env: + expansion_list.append(arr._extract(out_s)) + else: + data_map[out_s] = arr.env[out_s] + + elif isinstance(varr, np.ndarray): + data_map[arr_name] = varr + + elif varr is None: + pass # self + + elif np.isscalar(varr): + data_map[arr_name] = varr + + elif isinstance(varr, cl.array.Array): + data_map[arr_name] = varr + + else: + raise RuntimeError("unrecogonized captured variable %s" % arr_name) + + if knl is None: + # gather all known shapes + + for s in arr.shape: + if s in data_map: + continue + if arr._get_value(s) is not None: + data_map[s] = arr._get_value(s) + + for arg, _ in arr.namespace.values(): + if arg is None: + continue # skip self + assert isinstance(arg, Array) + for s in arg.shape: + if s not in arr.env: + expansion_list.append(arr._extract(s)) + elif arg._get_value(s) is not None: + data_map[s] = arr.env[s] + + else: + # try to get as much extra data that loopy wants as possible + # also evaluates the shape expressions + for argname, arg in knl.arg_dict.items(): + if argname in data_map: + continue + if argname in arr.env: + data_map[argname] = arr.env[argname] + elif argname in arr.namespace: + # Value unknown, but has expression. + # Try to expand the expression on inputs before runtime + # (applies to scalar expression, like shape vars) + expansion_list.append(arr._extract(argname)) + + # Make sure not to include None's prior to evaluating expressions + data_map = self._purge_nones(data_map) + + for se in expansion_list: + try: + seval = evaluate(se.expr, data_map) + except UnknownVariableError: + if 0: + warnings.warn( + "cannot get value for %s prior to calling " + "the loopy kernel" % se.name) + continue + + if se.name in data_map: + assert (seval is None) or (seval == data_map[se.name]) + else: + data_map[se.name] = seval + + # Remove keyward arguments not needed by the kernel + rm_args = [] + if knl is not None: + for argname in data_map.keys(): + if argname not in knl.arg_dict: + rm_args.append(argname) + for rmarg in rm_args: + del data_map[rmarg] + + return data_map + + def _purge_nones(self, data_map): + """Purge None-valued entries from a data map. + """ + entries_to_purge = [] + for key, val in data_map.items(): + if val is None: + entries_to_purge.append(key) + for key in entries_to_purge: + data_map.pop(key) + return data_map + + def _make_assignments(self, arr, **kwargs): + """Make the list of loopy instructions for a lazy array. + + :param substitutions: a dictionary of inlined information. + :param is_temporary: whether the lhs is a temporary. + """ + assignments = [] + + # a counter for generating unique names in diamond cases + sfx = '_%d' % self._counter + self._counter += 1 + + rhs_expr = arr.expr + if arr.ndim == 0: + lhs_expr = var(arr.name)[var(arr.name + '_dummy_iname')] + else: + lhs_expr = var(arr.name)[arr.inames] + + if 'substitutions' in kwargs: + rhs_expr = substitute(rhs_expr, kwargs['substitutions']) + + tc = TemporaryCollector(arr) + temps = tc(rhs_expr) + + is_temporary = arr.namespace[arr.name][1].get('is_temporary', False) + is_cse = arr.namespace[arr.name][1].get('is_cse', False) + has_barrier = arr.namespace[arr.name][1].get('has_barrier', False) + + if arr.name in temps: + # ignore self when dealing with temps of a temp + temps.remove(arr.name) + + if is_temporary: + assignments.append( + lp.Assignment( + lhs_expr, rhs_expr, + id=arr.name + sfx, + )) + if is_cse: + raise NotImplementedError() + if has_barrier: + raise NotImplementedError() + else: + assignments.append( + lp.Assignment(lhs_expr, rhs_expr,)) + + # recursively add assignments for all temps + if len(temps) > 0: + for tmp_name in temps: + tmp = arr._extract(tmp_name) + assignments.extend( + self._make_assignments(tmp, **kwargs)) + + return assignments + + def compile(self, *args, **kwargs): + # context: where the preconditions are checked + # data_map: passed to the loopy kernel + # scalar_map: inlined scalars + # + # context = data_map + scalar_map + context = dict() + + if len(args) > 1: + raise NotImplementedError() + arrs = list(args) + arr = arrs[0] + + shape_dtype = kwargs.pop('shape_dtype', np.int32) + + if arr.domain_expr is None: + loop_domain = arr._generate_index_domains() + else: + raise NotImplementedError() + + # collect kernel args + kernel_data = [] + if arr.ndim == 0: + kernel_data.append( + lp.GlobalArg( + arr.name, shape=(1, ), dtype=arr.dtype)) + else: + kernel_data.append( + lp.GlobalArg( + arr.name, shape=arr._shape_str, dtype=arr.dtype)) + + # evaluate scalar expressions + # e.g.: size of arange() + for aname, (a, atag) in arr.namespace.items(): + if a is None: + continue + if a.ndim == 0 and arr.env.get(aname, None) is None: + try: + env = arr.env.copy() + val = evaluate(a.expr, env) + if a.dtype is not None: + val = a.dtype(val) + arr.env[aname] = val + + except UnknownVariableError as e: + print(e) + pass + + inline_scalars = True + if inline_scalars: + scalar_map = dict() + assignments = self._make_assignments(arr, substitutions=scalar_map) + else: + assignments = self._make_assignments(arr) + + # add more kernel args (reduces a lot of guessing time) + # e.g., link array shape vars with the arrays, dtypes, etc. + # + # (the trick is to add no more than actually needed) + # FIXME: let loopy.ArgumentGuesser to find out what are needed. + for arr_name, (varr, tag) in arr.namespace.items(): + if arr_name == arr.name: + continue # arr (output) handled above + if isinstance(varr, Array) and varr.ndim > 0: + # for now, we only specify inputs + + if varr.dtype is None: + dtype = np.int32 if varr.is_integral else np.float64 + else: + dtype = varr.dtype + + if tag.get('is_argument', False): + kernel_data.append( + lp.GlobalArg( + arr_name, + dtype=dtype, + shape=varr._shape_str)) + elif tag.get('is_temporary', False): + kernel_data.append( + lp.TemporaryVariable( + arr_name, + dtype=dtype, + shape=varr._shape_str)) + + elif isinstance(varr, Array) and varr.ndim == 0: + if varr.dtype is None: + dtype = np.int32 if varr.is_integral else np.float64 + else: + dtype = varr.dtype + if tag.get('is_argument', False): + if arr_name in scalar_map: + pass + else: + kernel_data.append( + lp.ValueArg(arr_name, dtype=dtype)) + else: + pass + # warnings.warn( + # "cannot find shape information of %s" % arr_name) + + kernel_data.append('...') + + knl = lp.make_kernel( + loop_domain, assignments, + kernel_data, + lang_version=(2018, 2)) + + knl = lp.set_options(knl, return_dict=True) + + if inline_scalars: + for vn, (vv, vtag) in arr.namespace.items(): + if vv is None: + continue + if vv.ndim == 0: + scalar_map[vn] = arr._get_value(vn) + scalar_map = self._purge_nones(scalar_map) + context.update(scalar_map) + knl = lp.fix_parameters(knl, **scalar_map) + + # help with dtype inference + extra_dtype_dict = {} + for argname, arg in knl.arg_dict.items(): + if arg.dtype is None: + # array args + sym_arg = arr._extract(argname) + if sym_arg is not None: + extra_dtype_dict[argname] = sym_arg.dtype + + # shape args + if argname in arr.shape: + extra_dtype_dict[argname] = shape_dtype + for arr_arg in arr.bound_arguments.values(): + if arr_arg is None: + continue + if argname in arr_arg.shape: + extra_dtype_dict[argname] = shape_dtype + + knl = lp.add_and_infer_dtypes(knl, extra_dtype_dict) + + return knl, context + +# }}} End compiler diff --git a/lappy/eval/execution.py b/lappy/eval/execution.py new file mode 100644 index 0000000000000000000000000000000000000000..66f444f64e2ff55d95090bd3f960fb8c8b849265 --- /dev/null +++ b/lappy/eval/execution.py @@ -0,0 +1,36 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2020 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from lappy.eval import Closure + + +class Executor(object): + def __init__(self, closure): + assert isinstance(closure, Closure) + self.closure = closure + + def __call__(self, queue): + evt, lp_res = self.closure.kernel( + queue, **self.closure.data_map) + return lp_res diff --git a/lappy/exceptions.py b/lappy/exceptions.py new file mode 100644 index 0000000000000000000000000000000000000000..9829d0638881fdc29b9a431040ce1e12b9e821c6 --- /dev/null +++ b/lappy/exceptions.py @@ -0,0 +1,46 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2020 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +# {{{ errors + + +class AxisError(IndexError, ValueError): + # invalid axis passed + pass + +# }}} End errors + +# {{{ warnings + + +class ComplexWarning(RuntimeWarning): + # casting complex dtype to real + pass + + +class RankWarning(UserWarning): + # rank deficient matrix encountered when expecting otherwise + pass + +# }}} End warnings diff --git a/lappy/lib/__init__.py b/lappy/lib/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..6705d8fad967384b1b6d2cf9a33e333ad7542d22 --- /dev/null +++ b/lappy/lib/__init__.py @@ -0,0 +1,23 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2020 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" diff --git a/lappy/version.py b/lappy/version.py deleted file mode 100644 index 6ba2863610f2fe534ddbfabbbd4a41673be4294c..0000000000000000000000000000000000000000 --- a/lappy/version.py +++ /dev/null @@ -1,3 +0,0 @@ -VERSION = (2019, 1) -VERSION_STATUS = "" -VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 0000000000000000000000000000000000000000..c24fe5bb9e65c74aec64d71269f3617ed1bb776b --- /dev/null +++ b/pytest.ini @@ -0,0 +1,3 @@ +[pytest] +filterwarnings = + ignore::DeprecationWarning diff --git a/setup.py b/setup.py index 3fa8245a4fb482794425dad34168f01b73acf46c..f9171c05ee0811bfa9b9fdc01d59c1b9f56aef57 100644 --- a/setup.py +++ b/setup.py @@ -1,19 +1,58 @@ #! /usr/bin/env python # -*- coding: utf-8 -*- +import os +import subprocess from setuptools import setup, find_packages -ver_dic = {} -version_file = open("lappy/version.py") -try: - version_file_contents = version_file.read() -finally: - version_file.close() +YEAR = 2020 +MONTH = 1 +STATUS = 'dev' -exec(compile(version_file_contents, "lappy/version.py", 'exec'), ver_dic) + +def get_git_version(): + if not os.path.exists('.git'): + return None + + git_rev = None + try: + git_rev = subprocess.check_output( + ['git', 'rev-parse', 'HEAD'], stderr=subprocess.STDOUT + ).strip().decode('ascii') + except (subprocess.SubprocessError, OSError): + pass + + if git_rev is None: + git_rev = "Unknown" + + return git_rev + + +VERSION = (YEAR, MONTH) +VERSION_TEXT = '%d.%d%s' % (YEAR, MONTH, STATUS) +GIT_REVISION = get_git_version() + +VERSION_FILE = "lappy/version.py" +VER_DICT = { + 'version': VERSION, + 'version_status': STATUS, + 'version_text': VERSION_TEXT, + 'git_revision': GIT_REVISION} + +# rewrite version file every time +with open(VERSION_FILE, 'w') as a: + a.write("""################################################## +# THIS FILE IS GENERATED BY LAPPY SETUP.PY +################################################## + +VERSION = '%(version)s' +VERSION_STATUS = '%(version_status)s' +VERSION_TEXT = '%(version_text)s' +GIT_REVISION = '%(git_revision)s' +""" % VER_DICT) setup(name="lappy", - version="2017.1", + version=VER_DICT['version_text'], description="Lazy array programming in Python", long_description=open("README.rst").read(), classifiers=[ @@ -44,6 +83,11 @@ setup(name="lappy", "numpy", ], + extras_require={ + "test": ["pytest", ], + "doc": ["sphinx", "sphinx_rtd_theme"], + }, + author="Sophia Lin, Andreas Kloeckner, Xiaoyu Wei", url="http://pypi.python.org/pypi/lappy", author_email="svlin2@illinois.edu", diff --git a/test/test_array_api.py b/test/test_api_array.py similarity index 90% rename from test/test_array_api.py rename to test/test_api_array.py index 4bc703f103c62fe2b094cfc1a48835f55299cf7f..34adb410ed2301fed85e952f4312aa34f5c82fcb 100644 --- a/test/test_array_api.py +++ b/test/test_api_array.py @@ -23,7 +23,7 @@ THE SOFTWARE. """ import numpy as np -from lappy import Array, to_lappy_array +from lappy.core import Array, to_lappy_array import pytest @@ -36,7 +36,7 @@ import pytest def test_construction_from_numpy_array(nparr): arr = to_lappy_array(nparr) assert isinstance(arr, Array) - assert arr.shape == nparr.shape + assert arr._get_shape_vals_or_strs() == nparr.shape assert arr.dtype == nparr.dtype @@ -62,14 +62,14 @@ def test_construction_abstract(name, shape, dtype): arr = Array(name=name, shape=shape, dtype=dtype) assert arr.name == name assert arr.ndim == len(shape) - assert arr.shape == normalized_shape + assert arr._get_shape_vals_or_strs() == normalized_shape assert arr.dtype == dtype # name and shape are positional args arr = Array(name, shape, dtype=dtype) assert arr.name == name assert arr.ndim == len(shape) - assert arr.shape == normalized_shape + assert arr._get_shape_vals_or_strs() == normalized_shape assert arr.dtype == dtype # shape given by a string @@ -77,5 +77,5 @@ def test_construction_abstract(name, shape, dtype): arr = Array(name, shape_str, dtype=dtype) assert arr.name == name assert arr.ndim == len(shape) - assert arr.shape == normalized_shape + assert arr._get_shape_vals_or_strs() == normalized_shape assert arr.dtype == dtype diff --git a/test/test_broadcast.py b/test/test_broadcast.py index 92ebe9040f3bb7536f06fa8478e3b04fc9346058..1c62aca56f21c9092f7456ad472f6a5540a85541 100644 --- a/test/test_broadcast.py +++ b/test/test_broadcast.py @@ -25,8 +25,9 @@ THE SOFTWARE. import pytest import numpy as np from pymbolic import evaluate -from lappy.core.array import Array +from lappy import ndarray from lappy.core.broadcast import broadcast +from lappy.eval import Compiler @pytest.mark.parametrize('shapes', [ @@ -44,11 +45,11 @@ def test_broadcast_rules(shapes): nparrs = [np.zeros(s) for s in shapes] numpy_res = np.broadcast(*nparrs) - laarrs = [Array(shape=s) for s in shapes] + laarrs = [ndarray(shape=s) for s in shapes] lappy_res = broadcast(*laarrs) assert numpy_res.ndim == lappy_res.ndim - assert numpy_res.shape == lappy_res.shape + assert numpy_res.shape == lappy_res.shape_val_or_str assert numpy_res.size == lappy_res.size @@ -70,7 +71,7 @@ def test_symbolic_broadcast_rules(shapes): nparrs = [np.zeros(s) for s in shapes] numpy_res = np.broadcast(*nparrs) - laarrs = [Array(ndim=nd) for nd in ndims] + laarrs = [ndarray(ndim=nd) for nd in ndims] lappy_res = broadcast(*laarrs) assert numpy_res.ndim == lappy_res.ndim @@ -78,13 +79,21 @@ def test_symbolic_broadcast_rules(shapes): # inform lappy with actual shapes to evaluate exprs for i in range(len(laarrs)): laarrs[i] = laarrs[i].with_shape_data(shapes[i]) - arr_contexts = [arr.get_data_mapping() for arr in laarrs] + + compiler = Compiler() + closures = [compiler(arr) for arr in laarrs] + arr_contexts = [clos.data_map for clos in closures] broadcast_ctx = {} for ctx in arr_contexts: broadcast_ctx.update(ctx) + for arr in laarrs: + print(arr.env) + # shape and size are now expressions # evaluate those concretely with a sum + for lar in laarrs: + broadcast_ctx.update(lar.env) # check shape lappy_res_shape = [ diff --git a/test/test_indexing.py b/test/test_indexing.py new file mode 100644 index 0000000000000000000000000000000000000000..6cdd8b4da22f610a64a1e7e4c9d03727bdbedc5f --- /dev/null +++ b/test/test_indexing.py @@ -0,0 +1,45 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2020 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import lappy as la +import numpy as np +from lappy.core.basic_ops import SubArrayFactory + + +def test_subarr_out_ndim(): + arr = la.ndarray('A', ndim=2) + subfac = SubArrayFactory( + arr, + (slice(None, None, None), [[0, 1], [1, 1]], np.newaxis)) + print(subfac.selectors) + print(subfac.ndim) + print(arr.inames) + print(arr.shape) + + +def test_slice_completion(): + arr = la.ndarray('A', ndim=2) + subfac = SubArrayFactory(arr) + sl = subfac._complete_slice(slice(None, -10, -1), 0) + assert sl.step == -1 diff --git a/test/test_numpy_consistency.py b/test/test_numpy_consistency.py new file mode 100644 index 0000000000000000000000000000000000000000..ddec2f5779282d008ff2aedb459dd84a1e6c2170 --- /dev/null +++ b/test/test_numpy_consistency.py @@ -0,0 +1,123 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy # noqa: F401 +import lappy # noqa: F401 + +import pytest +from warnings import warn + +# part of the numpy api is not meant to be included. +NOSUP = [ + '__array_finalize__', '__array_prepare__', + '__array_wrap__', + '__getstate__', '__setstate__', '__array_priority__', + '__array_struct__', '__long__', '__setmask__', + '_baseclass', '_comparison', '_data', '_defaulthardmask', + '_defaultmask', '_delegate_binop', '_get_data', '_get_flat', + '_get_mask', '_get_recordmask', '_insert_masked_print', '_print_width', + '_print_width_1d', '_set_flat', '_set_mask', '_set_recordmask', + '_update_from', + + # build system + '__config__', '__NUMPY_SETUP__', '_NoValue', 'MachAr', 'WRAP', + 'Tester', 'ModuleDeprecationWarning', '_add_newdoc_ufunc', '_arg', + '_distributor_init', '_pytesttester' + + # not well-defined for lappy arrays, since there is no direct + # correspondence with data in memory. + 'iscontiguous', 'byteswap', 'newbyteorder', 'getfield', 'setfield', + 'strides', 'tobytes', 'tofile', 'tolist', 'tostring', 'view', + 'ALLOW_THREADS', 'BUFSIZE', 'CLIP', 'ERR_CALL', 'ERR_DEFAULT', + 'ERR_IGNORE', 'ERR_LOG', 'ERR_PRINT', 'ERR_RAISE', 'ERR_WARN', + 'FLOATING_POINT_SUPPORT', 'FPE_DIVIDEBYZERO', 'FPE_INVALID', + 'FPE_OVERFLOW', 'FPE_UNDERFLOW', 'False_', 'True_', 'RAISE', + 'MAXDIMS', 'MAY_SHARE_BOUNDS', 'MAY_SHARE_EXACT', + 'Inf', 'Infinity', 'NAN', 'NINF', 'NZERO', 'NaN', 'PINF', 'PZERO', + 'SHIFT_DIVIDEBYZERO', 'SHIFT_INVALID', 'SHIFT_OVERFLOW', 'SHIFT_UNDERFLOW', + 'UFUNC_BUFSIZE_DEFAULT', 'UFUNC_PYVALS_NAME', '_UFUNC_API', + + # not planned to support + 'DataSource', '_globals', + + # (no masked array support for now) + 'compress', 'compressed', 'count', 'base', 'baseclass', + 'choose', 'ctypes', 'data', 'fill_value', 'filled', + 'get_fill_value', 'get_imag', 'get_real', + 'harden_mask', 'hardmask', 'ids', 'mini', 'mask', 'recordmask', + 'set_fill_value', 'sharedmask', 'shrink_mask', 'soften_mask', + 'toflex', 'torecords', 'unshare_mask', + + # (no matrixlib support for now) + '_mat', + + # attributes set during construction + 'dtype', 'shape', 'flags', 'setflags', + ] + + +def warn_diff(npi, lpi, name, assert_extra_hidden=False): + missing = [] + extra = [] + + for fn in npi: + if fn not in lpi: + missing.append(fn) + + missing = [m for m in missing if m not in NOSUP] + + for fn in lpi: + if fn not in npi: + extra.append(fn) + + if len(missing) > 0: + warn("\n>>>> Missing API in %s (%d):\n %s" + % (name, len(missing), ', '.join(missing))) + if len(extra) > 0: + warn("\n>>>> Extra API in %s (%d):\n %s" + % (name, len(extra), ', '.join(extra))) + + if assert_extra_hidden: + assert all(ext.startswith('_') for ext in extra) + + return len(missing) + len(extra) + + +@pytest.mark.parametrize('path_lappy, path_numpy, name, tol', [ + ('lappy', 'numpy', 'lappy', 1000), + ('lappy.ndarray', 'numpy.ndarray, numpy.ma.MaskedArray', 'ndarray', 1000), + ]) +def test_public_interface(path_lappy, path_numpy, name, tol): + lpi = [] + npi = [] + for pl in path_lappy.split(', '): + exec('lpi += dir(%s)' % pl) + for pn in path_numpy.split(', '): + exec('npi += dir(%s)' % pn) + + lpi = sorted(list(set(lpi))) + npi = sorted(list(set(npi))) + + distance = warn_diff(npi, lpi, name) + assert distance < tol diff --git a/test/test_reshape.py b/test/test_reshape.py index 66eb8bbfb6eeca8e477965a234529a033dac85ae..d117e8270f9ae0fdf629743966b6a016042efada 100644 --- a/test/test_reshape.py +++ b/test/test_reshape.py @@ -23,6 +23,7 @@ THE SOFTWARE. """ import pytest +import pyopencl as cl import lappy as la import numpy as np @@ -49,21 +50,18 @@ from pyopencl.tools import ( # noqa ]) @pytest.mark.parametrize('order', ['C', 'F']) def test_reshape(ctx_factory, in_shape, out_shape, order, dtype=np.float64): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) ndim = len(in_shape) - env = la.core.array.default_env() - env['cl_ctx'] = ctx_factory() - - A = la.Array('A', ndim=ndim, env=env) # noqa: N806 + A = la.ndarray('A', ndim=ndim) # noqa: N806 B = A.reshape(out_shape, order=order, name='B') # noqa: N806 data = np.random.rand(*in_shape).astype(dtype) C = B.with_data(A=data).with_name('C') # noqa: N806 - print(C.preconditions) - - lappy_val = C.eval() + lappy_val = C.eval(queue) numpy_val = data.reshape(out_shape, order=order) assert np.allclose(lappy_val, numpy_val) diff --git a/test/test_transpose.py b/test/test_transpose.py index 3ea5516130752a3ce7e56f6b887821de161d2500..8ba0b4c10e7d478479b7a780023eb62ec54ce4aa 100644 --- a/test/test_transpose.py +++ b/test/test_transpose.py @@ -23,6 +23,7 @@ THE SOFTWARE. """ import pytest +import pyopencl as cl import lappy as la import numpy as np @@ -37,21 +38,18 @@ from pyopencl.tools import ( # noqa ((2, 3, 4, 6), (0, 2, 1, 3)), ]) def test_transpose(ctx_factory, test_shape, axes, dtype=np.float64): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) ndim = len(test_shape) - env = la.core.array.default_env() - env['cl_ctx'] = ctx_factory() - - A = la.Array('A', ndim=ndim, env=env) # noqa: N806 + A = la.ndarray('A', ndim=ndim) # noqa: N806 B = A.transpose(axes=axes, name='B') # noqa: N806 data = np.random.rand(*test_shape).astype(dtype) C = B.with_data(A=data).with_name('C') # noqa: N806 - print(repr(C)) - print(C.inames) - lappy_val = C.eval() + lappy_val = C.eval(queue) numpy_val = data.transpose(axes) assert np.allclose(lappy_val, numpy_val) diff --git a/test/test_ufunc.py b/test/test_ufunc.py index 296cc8b885c85aa421ed67175a8ae7b4f75ded79..8edb082e49422eb470e9ec172323bdbbf085a925 100644 --- a/test/test_ufunc.py +++ b/test/test_ufunc.py @@ -24,6 +24,7 @@ THE SOFTWARE. import numpy as np import lappy as la +import pyopencl as cl import lappy.core.ufuncs as math import pytest @@ -38,12 +39,12 @@ from pyopencl.tools import ( # noqa ('exp', [(2, 4, 8), (1, 2, 3), ], np.complex128), ]) def test_unary_ufunc(ctx_factory, ufunc, shapes, dtype): - env = la.core.array.default_env() - env['cl_ctx'] = ctx_factory() + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) ndim = len(shapes[0]) # symbolic code is reused for different shapes (same ndim) - mat = la.Array('A', ndim=ndim, dtype=dtype, env=env) + mat = la.ndarray('A', ndim=ndim, dtype=dtype) fmat = getattr(math, ufunc)(mat).with_name('B') for shape in shapes: @@ -52,7 +53,7 @@ def test_unary_ufunc(ctx_factory, ufunc, shapes, dtype): bound_mat = fmat.with_data(A=mat_data).with_name('C') # eval - lappy_res = bound_mat.eval() + lappy_res = bound_mat.eval(queue) numpy_res = getattr(np, ufunc)(mat_data) # compare with numpy @@ -65,14 +66,14 @@ def test_unary_ufunc(ctx_factory, ufunc, shapes, dtype): ('multiply', [(2, 4, 8), (1, 2, 3), ], np.complex128), ]) def test_binary_ufunc(ctx_factory, ufunc, shapes, dtype): - env = la.core.array.default_env() - env['cl_ctx'] = ctx_factory() + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) ndim = len(shapes[0]) sym_shape = tuple('s%d' % i for i in range(ndim)) - mat_a = la.Array('A', shape=sym_shape, dtype=dtype, env=env) - mat_b = la.Array('B', shape=sym_shape, dtype=dtype, env=env) + mat_a = la.ndarray('A', shape=sym_shape, dtype=dtype) + mat_b = la.ndarray('B', shape=sym_shape, dtype=dtype) fmat = getattr(math, ufunc)(mat_a, mat_b).with_name('C') @@ -87,7 +88,7 @@ def test_binary_ufunc(ctx_factory, ufunc, shapes, dtype): ).with_name('C') # eval - lappy_res = bound_mat.eval() + lappy_res = bound_mat.eval(queue) numpy_res = getattr(np, ufunc)(mat_a_data, mat_b_data) # compare with numpy @@ -106,16 +107,16 @@ def test_binary_ufunc(ctx_factory, ufunc, shapes, dtype): @pytest.mark.parametrize('ufunc', ['add', 'subtract', 'multiply', 'divide']) @pytest.mark.parametrize('dtype', [np.float64, ]) def test_binary_ufunc_with_broadcast(ctx_factory, ufunc, shapes, dtype): - env = la.core.array.default_env() - env['cl_ctx'] = ctx_factory() + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) ndim_a = len(shapes[0]) ndim_b = len(shapes[1]) sym_shape_a = tuple('s_a_%d' % i for i in range(ndim_a)) sym_shape_b = tuple('s_b_%d' % i for i in range(ndim_b)) - mat_a = la.Array('A', shape=sym_shape_a, dtype=dtype, env=env) - mat_b = la.Array('B', shape=sym_shape_b, dtype=dtype, env=env) + mat_a = la.ndarray('A', shape=sym_shape_a, dtype=dtype) + mat_b = la.ndarray('B', shape=sym_shape_b, dtype=dtype) fmat = getattr(math, ufunc)(mat_a, mat_b).with_name('C') @@ -129,7 +130,7 @@ def test_binary_ufunc_with_broadcast(ctx_factory, ufunc, shapes, dtype): ).with_name('C') # eval - lappy_res = bound_mat.eval() + lappy_res = bound_mat.eval(queue) numpy_res = getattr(np, ufunc)(mat_a_data, mat_b_data) # compare with numpy