BUG: .dot() not producing the same as sparse.basic.dot()

mattijsdp commented 1 year ago

Describe the issue:

When I create a sparse matrix and multiply this with a vector using the class' method I get a NotImplementedError whereas when I use pytensor.sparse.basic.dot() it does work.

Reproducable code example:

import pymc as pm
import numpy as np
from pytensor import sparse
y = np.random.randn(3)

data = np.asarray([7, 8, 9])
indices = np.asarray([0, 1, 2])
indptr = np.asarray([0, 2, 3, 3])
X_sparse = sparse.CSC(data=data, indices=indices, indptr=indptr, shape=(3, 3))

with pm.Model() as py_lr2:
    intercept = pm.Normal('c', 0, 100)
    betas = pm.Normal('betas', 0, 10, shape=(3))
    sd = pm.HalfNormal('sigma', 5.)

    # Option 1: this compiles
    # mu = sparse.basic.dot(X_sparse, betas)

    # Option 2: this raises a NotImplementedError
    mu = X_sparse.dot(betas)

    obs = pm.Normal('y', mu=mu + intercept, sigma=sd, observed=y)

Error message:

NotImplementedError                       Traceback (most recent call last)
Cell In[2], line 20
     14 sd = pm.HalfNormal('sigma', 5.)
     16 # Option 1: this compiles
     17 # mu = sparse.basic.dot(X_sparse, betas)
     18 
     19 # Option 2: this raises a NotImplementedError
---> 20 mu = X_sparse.dot(betas)
     22 obs = pm.Normal('y', mu=mu + intercept, sigma=sd, observed=y)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:375, in _sparse_py_operators.__dot__(left, right)
    374 def __dot__(left, right):
--> 375     return structured_dot(left, right)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:3557, in structured_dot(x, y)
   3554     raise TypeError("structured_dot requires at least one sparse argument")
   3556 if x_is_sparse_variable:
-> 3557     return _structured_dot(x, y)
   3558 else:
   3559     assert y_is_sparse_variable

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/graph/op.py:295, in Op.__call__(self, *inputs, **kwargs)
    253 r"""Construct an `Apply` node using :meth:`Op.make_node` and return its outputs.
    254 
    255 This method is just a wrapper around :meth:`Op.make_node`.
   (...)
    292 
    293 """
    294 return_list = kwargs.pop("return_list", False)
--> 295 node = self.make_node(*inputs, **kwargs)
    297 if config.compute_test_value != "off":
    298     compute_test_value(node)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:3437, in StructuredDot.make_node(self, a, b)
   3435 dtype_out = aes.upcast(a.type.dtype, b.type.dtype)
   3436 if b.type.ndim != 2:
-> 3437     raise NotImplementedError("non-matrix b")
   3439 if _is_sparse_variable(b):
   3440     return Apply(self, [a, b], [SparseTensorType(a.type.format, dtype_out)()])

NotImplementedError: non-matrix b

PyTensor version information:

Operating system: linux Conda installation (channel conda-forge):

Python version: 3.10.9
PyTensor version: 2.12.1
Pymc version: 5.4.0
Numpy version: 1.24.3
floatX ({'float32', 'float16', 'float64'}) Doc: Default floating-point precision for python casts.

Note: float16 support is experimental, use at your own risk. Value: float64

warn_float64 ({'pdb', 'raise', 'ignore', 'warn'}) Doc: Do an action when a tensor variable with float64 dtype is created. Value: ignore

pickle_test_value (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041840>>) Doc: Dump test values while pickling model. If True, test values will be dumped with model. Value: True

cast_policy ({'numpy+floatX', 'custom'}) Doc: Rules for implicit type casting Value: custom

deterministic ({'default', 'more'}) Doc: If more, sometimes we will select some implementation that are more deterministic, but slower. Also see the dnn.conv.algo* flags to cover more cases. Value: default

device (cpu) Doc: Default device for computations. only cpu is supported for now Value: cpu

force_device (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340418d0>>) Doc: Raise an error if we can't use the specified device Value: False

conv__assert_shape (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041900>>) Doc: If True, AbstractConv* ops will verify that user-provided shapes match the runtime shapes (debugging option, may slow down compilation) Value: False

print_global_stats (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041930>>) Doc: Print some global statistics (time spent) at the end Value: False

assert_no_cpu_op ({'pdb', 'raise', 'ignore', 'warn'}) Doc: Raise an error/warning if there is a CPU op in the computational graph. Value: ignore

unpickle_function (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041990>>) Doc: Replace unpickled PyTensor functions with None. This is useful to unpickle old graphs that pickled them when it shouldn't Value: True

<pytensor.configparser.ConfigParam object at 0x7ff2340419c0> Doc: Default compilation mode Value: Mode

cxx (<class 'str'>) Doc: The C++ compiler to use. Currently only g++ is supported, but supporting additional compilers should not be too difficult. If it is empty, no C++ code is compiled. Value: /usr/bin/g++

linker ({'c|py', 'vm_nogc', 'cvm', 'c', 'cvm_nogc', 'vm', 'py', 'c|py_nogc'}) Doc: Default linker used if the pytensor flags mode is Mode Value: cvm

allow_gc (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041c00>>) Doc: Do we default to delete intermediate results during PyTensor function calls? Doing so lowers the memory requirement, but asks that we reallocate memory at the next function call. This is implemented for the default linker, but may not work for all linkers. Value: True

optimizer ({'o1', 'o3', 'merge', 'o2', 'o4', 'fast_run', 'None', 'fast_compile', 'unsafe'}) Doc: Default optimizer. If not None, will use this optimizer with the Mode Value: o4

optimizer_verbose (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041a20>>) Doc: If True, we print all optimization being applied Value: False

on_opt_error ({'pdb', 'raise', 'ignore', 'warn'}) Doc: What to do when an optimization crashes: warn and skip it, raise the exception, or fall into the pdb debugger. Value: warn

nocleanup (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041c30>>) Doc: Suppress the deletion of code files that did not compile cleanly Value: False

on_unused_input ({'warn', 'ignore', 'raise'}) Doc: What to do if a variable in the 'inputs' list of pytensor.function() is not used in the graph. Value: raise

gcc__cxxflags (<class 'str'>) Doc: Extra compiler flags for gcc Value: -Wno-c++11-narrowing -fno-exceptions -fno-unwind-tables -fno-asynchronous-unwind-tables

cmodule__warn_no_version (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041ab0>>) Doc: If True, will print a warning when compiling one or more Op with C code that can't be cached because there is no c_code_cache_version() function associated to at least one of those Ops. Value: False

cmodule__remove_gxx_opt (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041b10>>) Doc: If True, will remove the -O* parameter passed to g++.This is useful to debug in gdb modules compiled by PyTensor.The parameter -g is passed by default to g++ Value: False

cmodule__compilation_warning (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041b70>>) Doc: If True, will print compilation warnings. Value: False

cmodule__preload_cache (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041ae0>>) Doc: If set to True, will preload the C module cache at import time Value: False

cmodule__age_thresh_use (<class 'int'>) Doc: In seconds. The time after which PyTensor won't reuse a compile c module. Value: 2073600

cmodule__debug (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041cc0>>) Doc: If True, define a DEBUG macro (if not exists) for any compiled C code. Value: False

compile__wait (<class 'int'>) Doc: Time to wait before retrying to acquire the compile lock. Value: 5

compile__timeout (<class 'int'>) Doc: In seconds, time that a process will wait before deciding to override an existing lock. An override only happens when the existing lock is held by the same owner and has not been 'refreshed' by this owner for more than this period. Refreshes are done every half timeout period for running processes. Value: 120

ctc__root (<class 'str'>) Doc: Directory which contains the root of Baidu CTC library. It is assumed that the compiled library is either inside the build, lib or lib64 subdirectory, and the header inside the include directory. Value:

tensor__cmp_sloppy (<class 'int'>) Doc: Relax pytensor.tensor.math._allclose (0) not at all, (1) a bit, (2) more Value: 0

tensor__local_elemwise_fusion (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041db0>>) Doc: Enable or not in fast_run mode(fast_run optimization) the elemwise fusion optimization Value: True

lib__amblibm (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041de0>>) Doc: Use amd's amdlibm numerical library Value: False

tensor__insert_inplace_optimizer_validate_nb (<class 'int'>) Doc: -1: auto, if graph have less then 500 nodes 1, else 10 Value: -1

traceback__limit (<class 'int'>) Doc: The number of stack to trace. -1 mean all. Value: 8

traceback__compile_limit (<class 'int'>) Doc: The number of stack to trace to keep during compilation. -1 mean all. If greater then 0, will also make us save PyTensor internal stack trace. Value: 0

experimental__local_alloc_elemwise (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041ea0>>) Doc: DEPRECATED: If True, enable the experimental optimization local_alloc_elemwise. Generates error if not True. Use optimizer_excluding=local_alloc_elemwise to disable. Value: True

experimental__local_alloc_elemwise_assert (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041ed0>>) Doc: When the local_alloc_elemwise is applied, add an assert to highlight shape errors. Value: True

warn__ignore_bug_before ({'1.0.4', '0.8.1', '0.4.1', '0.3', '0.8', '1.0.2', '0.5', 'None', '1.0.5', '0.8.2', '1.0.3', '0.7', '0.4', 'all', '1.0.1', '0.6', '0.10', '1.0', '0.9'}) Doc: If 'None', we warn about all PyTensor bugs found by default. If 'all', we don't warn about PyTensor bugs found by default. If a version, we print only the warnings relative to PyTensor bugs found after that version. Warning for specific bugs can be configured with specific [warn] flags. Value: 0.9

exception_verbosity ({'high', 'low'}) Doc: If 'low', the text of exceptions will generally refer to apply nodes with short names such as Elemwise{add_no_inplace}. If 'high', some exceptions will also refer to apply nodes with long descriptions like: A. Elemwise{add_no_inplace} B. log_likelihood_v_given_h C. log_likelihood_h Value: low

print_test_value (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041f60>>) Doc: If 'True', the eval of an PyTensor variable will return its test_value when this is available. This has the practical consequence that, e.g., in debugging my_var will print the same as my_var.tag.test_value when a test value is defined. Value: False

compute_test_value ({'pdb', 'ignore', 'raise', 'off', 'warn'}) Doc: If 'True', PyTensor will run each op at graph build time, using Constants, SharedVariables and the tag 'test_value' as inputs to the function. This helps the user track down problems in the graph before it gets optimized. Value: off

compute_test_value_opt ({'pdb', 'ignore', 'raise', 'off', 'warn'}) Doc: For debugging PyTensor optimization only. Same as compute_test_value, but is used during PyTensor optimization Value: off

check_input (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234041ff0>>) Doc: Specify if types should check their input in their C code. It can be used to speed up compilation, reduce overhead (particularly for scalars) and reduce the number of generated C files. Value: True

NanGuardMode__nan_is_error (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042020>>) Doc: Default value for nan_is_error Value: True

NanGuardMode__inf_is_error (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042050>>) Doc: Default value for inf_is_error Value: True

NanGuardMode__big_is_error (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042080>>) Doc: Default value for big_is_error Value: True

NanGuardMode__action ({'pdb', 'warn', 'raise'}) Doc: What NanGuardMode does when it finds a problem Value: raise

DebugMode__patience (<class 'int'>) Doc: Optimize graph this many times to detect inconsistency Value: 10

DebugMode__check_c (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042110>>) Doc: Run C implementations where possible Value: True

DebugMode__check_py (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042140>>) Doc: Run Python implementations where possible Value: True

DebugMode__check_finite (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042170>>) Doc: True -> complain about NaN/Inf results Value: True

DebugMode__check_strides (<class 'int'>) Doc: Check that Python- and C-produced ndarrays have same strides. On difference: (0) - ignore, (1) warn, or (2) raise error Value: 0

DebugMode__warn_input_not_reused (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340421d0>>) Doc: Generate a warning when destroy_map or view_map says that an op works inplace, but the op did not reuse the input for its output. Value: True

DebugMode__check_preallocated_output (<class 'str'>) Doc: Test thunks with pre-allocated memory as output storage. This is a list of strings separated by ":". Valid values are: "initial" (initial storage in storage map, happens with Scan),"previous" (previously-returned memory), "c_contiguous", "f_contiguous", "strided" (positive and negative strides), "wrong_size" (larger and smaller dimensions), and "ALL" (all of the above). Value:

DebugMode__check_preallocated_output_ndim (<class 'int'>) Doc: When testing with "strided" preallocated output memory, test all combinations of strides over that number of (inner-most) dimensions. You may want to reduce that number to reduce memory or time usage, but it is advised to keep a minimum of 2. Value: 4

profiling__time_thunks (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042260>>) Doc: Time individual thunks when profiling Value: True

profiling__n_apply (<class 'int'>) Doc: Number of Apply instances to print by default Value: 20

profiling__n_ops (<class 'int'>) Doc: Number of Ops to print by default Value: 20

profiling__output_line_width (<class 'int'>) Doc: Max line width for the profiling output Value: 512

profiling__min_memory_size (<class 'int'>) Doc: For the memory profile, do not print Apply nodes if the size of their outputs (in bytes) is lower than this threshold Value: 1024

profiling__min_peak_memory (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042350>>) Doc: The min peak memory usage of the order Value: False

profiling__destination (<class 'str'>) Doc: File destination of the profiling output Value: stderr

profiling__debugprint (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340423b0>>) Doc: Do a debugprint of the profiled functions Value: False

profiling__ignore_first_call (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340423e0>>) Doc: Do we ignore the first call of an PyTensor function. Value: False

on_shape_error ({'raise', 'warn'}) Doc: warn: print a warning and use the default value. raise: raise an error Value: warn

openmp (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042440>>) Doc: Allow (or not) parallel computation on the CPU with OpenMP. This is the default value used when creating an Op that supports OpenMP parallelization. It is preferable to define it via the PyTensor configuration file ~/.pytensorrc or with the environment variable PYTENSOR_FLAGS. Parallelization is only done for some operations that implement it, and even for operations that implement parallelism, each operation is free to respect this flag or not. You can control the number of threads used with the environment variable OMP_NUM_THREADS. If it is set to 1, we disable openmp in PyTensor by default. Value: False

openmp_elemwise_minsize (<class 'int'>) Doc: If OpenMP is enabled, this is the minimum size of vectors for which the openmp parallelization is enabled in element wise ops. Value: 200000

optimizer_excluding (<class 'str'>) Doc: When using the default mode, we will remove optimizer with these tags. Separate tags with ':'. Value:

optimizer_including (<class 'str'>) Doc: When using the default mode, we will add optimizer with these tags. Separate tags with ':'. Value:

optimizer_requiring (<class 'str'>) Doc: When using the default mode, we will require optimizer with these tags. Separate tags with ':'. Value:

optdb__position_cutoff (<class 'float'>) Doc: Where to stop earlier during optimization. It represent the position of the optimizer where to stop. Value: inf

optdb__max_use_ratio (<class 'float'>) Doc: A ratio that prevent infinite loop in EquilibriumGraphRewriter. Value: 8.0

cycle_detection ({'regular', 'fast'}) Doc: If cycle_detection is set to regular, most inplaces are allowed,but it is slower. If cycle_detection is set to faster, less inplacesare allowed, but it makes the compilation faster.The interaction of which one give the lower peak memory usage iscomplicated and not predictable, so if you are close to the peakmemory usage, triyng both could give you a small gain. Value: regular

check_stack_trace ({'warn', 'log', 'raise', 'off'}) Doc: A flag for checking the stack trace during the optimization process. default (off): does not check the stack trace of any optimization log: inserts a dummy stack trace that identifies the optimizationthat inserted the variable that had an empty stack trace.warn: prints a warning if a stack trace is missing and also a dummystack trace is inserted that indicates which optimization insertedthe variable that had an empty stack trace.raise: raises an exception if a stack trace is missing Value: off

metaopt__verbose (<class 'int'>) Doc: 0 for silent, 1 for only warnings, 2 for full output withtimings and selected implementation Value: 0

metaopt__optimizer_excluding (<class 'str'>) Doc: exclude optimizers with these tags. Separate tags with ':'. Value:

metaopt__optimizer_including (<class 'str'>) Doc: include optimizers with these tags. Separate tags with ':'. Value:

unittests__rseed (<class 'str'>) Doc: Seed to use for randomized unit tests. Special value 'random' means using a seed of None. Value: 666

warn__round (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340426b0>>) Doc: Warn when using tensor.round with the default mode. Round changed its default from half_away_from_zero to half_to_even to have the same default as NumPy. Value: False

profile (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340426e0>>) Doc: If VM should collect profile information Value: False

profile_optimizer (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042710>>) Doc: If VM should collect optimizer profile information Value: False

profile_memory (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042740>>) Doc: If VM should collect memory profile information and print it Value: False

<pytensor.configparser.ConfigParam object at 0x7ff234042770> Doc: Useful only for the VM Linkers. When lazy is None, auto detect if lazy evaluation is needed and use the appropriate version. If the C loop isn't being used and lazy is True, use the Stack VM; otherwise, use the Loop VM. Value: None

numba__vectorize_target ({'cpu', 'cuda', 'parallel'}) Doc: Default target for numba.vectorize. Value: cpu

numba__fastmath (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff2340427d0>>) Doc: If True, use Numba's fastmath mode. Value: True

numba__cache (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff234042800>>) Doc: If True, use Numba's file based caching. Value: True

compiledir_format (<class 'str'>) Doc: Format string for platform-dependent compiled module subdirectory (relative to base_compiledir). Available keys: device, gxx_version, hostname, numpy_version, platform, processor, pytensor_version, python_bitwidth, python_int_bitwidth, python_version, shortplatform. Defaults to compiledir%(short_platform)s-%(processor)s- %(python_version)s-%(pythonbitwidth)s. Value: compiledir%(short_platform)s-%(processor)s-%(python_version)s-%(python_bitwidth)s

<pytensor.configparser.ConfigParam object at 0x7ff2340428c0> Doc: platform-independent root directory for compiled modules Value: /home/jupyter-mattijs/.pytensor

<pytensor.configparser.ConfigParam object at 0x7ff234042830> Doc: platform-dependent cache directory for compiled modules Value: /home/jupyter-mattijs/.pytensor/compiledir_Linux-5.15--aws-x86_64-with-glibc2.31-x86_64-3.10.9-64

blas__ldflags (<class 'str'>) Doc: lib[s] to include for [Fortran] level-3 blas implementation Value: -L/home/jupyter-mattijs/.conda/envs/dwell-times/lib -lmkl_core -lmkl_intel_thread -lmkl_rt -Wl,-rpath,/home/jupyter-mattijs/.conda/envs/dwell-times/lib

blas__check_openmp (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff1ad60be80>>) Doc: Check for openmp library conflict. WARNING: Setting this to False leaves you open to wrong results in blas-related operations. Value: True

scan__allow_gc (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff1a0fad9c0>>) Doc: Allow/disallow gc inside of Scan (default: False) Value: False

scan__allow_output_prealloc (<bound method BoolParam._apply of <pytensor.configparser.BoolParam object at 0x7ff1aaec6980>>) Doc: Allow/disallow memory preallocation for outputs inside of scan (default: True) Value: True

Context for the issue:

I think it is simply counterintuitive and don't understand why .dot() raises an error. Admittedly, I am very much a beginner when it comes to PyTensor so I might be missing something.

ricardoV94 commented 1 year ago

It seems to try to dispatch to structure_dot, which must be matrix-matrix. We should fallback to dot when we have matrix-vector?

mattijsdp commented 1 year ago

Turns out that although using .dot() allows you to create the model, the resulting graph seems to have problems as well as it leads to the same error when sampling: i.e. running idata = pm.sample(model=py_lr2) leads to:


File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/tensor/math.py:2053, in dot(l, r)
   2052 try:
-> 2053     res = l.__dot__(r)
   2054     if res is NotImplemented:

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:375, in _sparse_py_operators.__dot__(left, right)
    374 def __dot__(left, right):
--> 375     return structured_dot(left, right)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:3557, in structured_dot(x, y)
   3556 if x_is_sparse_variable:
-> 3557     return _structured_dot(x, y)
   3558 else:

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/graph/op.py:295, in Op.__call__(self, *inputs, **kwargs)
    294 return_list = kwargs.pop("return_list", False)
--> 295 node = self.make_node(*inputs, **kwargs)
    297 if config.compute_test_value != "off":

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:3437, in StructuredDot.make_node(self, a, b)
   3436 if b.type.ndim != 2:
-> 3437     raise NotImplementedError("non-matrix b")
   3439 if _is_sparse_variable(b):

NotImplementedError: non-matrix b

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[68], line 1
----> 1 idata = pm.sample(model=py_lr2)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pymc/sampling/mcmc.py:634, in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
    631         auto_nuts_init = False
    633 initial_points = None
--> 634 step = assign_step_methods(model, step, methods=pm.STEP_METHODS, step_kwargs=kwargs)
    636 if nuts_sampler != "pymc":
    637     if not isinstance(step, NUTS):

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pymc/sampling/mcmc.py:200, in assign_step_methods(model, step, methods, step_kwargs)
    198 if has_gradient:
    199     try:
--> 200         tg.grad(model_logp, var)
    201     except (NotImplementedError, tg.NullTypeGradError):
    202         has_gradient = False

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/gradient.py:617, in grad(cost, wrt, consider_constant, disconnected_inputs, add_names, known_grads, return_disconnected, null_gradients)
    614     if hasattr(g.type, "dtype"):
    615         assert g.type.dtype in pytensor.tensor.type.float_dtypes
--> 617 _rval: Sequence[Variable] = _populate_grad_dict(
    618     var_to_app_to_idx, grad_dict, _wrt, cost_name
    619 )
    621 rval: MutableSequence[Optional[Variable]] = list(_rval)
    623 for i in range(len(_rval)):

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/gradient.py:1420, in _populate_grad_dict(var_to_app_to_idx, grad_dict, wrt, cost_name)
   1417     # end if cache miss
   1418     return grad_dict[var]
-> 1420 rval = [access_grad_cache(elem) for elem in wrt]
   1422 return rval

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/gradient.py:1420, in <listcomp>(.0)
   1417     # end if cache miss
   1418     return grad_dict[var]
-> 1420 rval = [access_grad_cache(elem) for elem in wrt]
   1422 return rval

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/gradient.py:1375, in _populate_grad_dict.<locals>.access_grad_cache(var)
   1373 for node in node_to_idx:
   1374     for idx in node_to_idx[node]:
-> 1375         term = access_term_cache(node)[idx]
   1377         if not isinstance(term, Variable):
   1378             raise TypeError(
   1379                 f"{node.op}.grad returned {type(term)}, expected"
   1380                 " Variable instance."
   1381             )

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/gradient.py:1205, in _populate_grad_dict.<locals>.access_term_cache(node)
   1197         if o_shape != g_shape:
   1198             raise ValueError(
   1199                 "Got a gradient of shape "
   1200                 + str(o_shape)
   1201                 + " on an output of shape "
   1202                 + str(g_shape)
   1203             )
-> 1205 input_grads = node.op.L_op(inputs, node.outputs, new_output_grads)
   1207 if input_grads is None:
   1208     raise TypeError(
   1209         f"{node.op}.grad returned NoneType, expected iterable."
   1210     )

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/graph/op.py:392, in Op.L_op(self, inputs, outputs, output_grads)
    365 def L_op(
    366     self,
    367     inputs: Sequence[Variable],
    368     outputs: Sequence[Variable],
    369     output_grads: Sequence[Variable],
    370 ) -> List[Variable]:
    371     r"""Construct a graph for the L-operator.
    372 
    373     The L-operator computes a row vector times the Jacobian.
   (...)
    390 
    391     """
--> 392     return self.grad(inputs, output_grads)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/sparse/basic.py:4025, in Dot.grad(self, inputs, gout)
   4023     rval.append(dot(gz, y.T))
   4024 if _is_dense_variable(x):
-> 4025     rval.append(at_dot(x.T, gz))
   4026 else:
   4027     rval.append(dot(x.T, gz))

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/tensor/math.py:2057, in dot(l, r)
   2055         raise NotImplementedError
   2056 except (NotImplementedError, AttributeError, TypeError):
-> 2057     res = r.__rdot__(l)
   2058     if res is NotImplemented:
   2059         raise NotImplementedError()

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/tensor/var.py:650, in _tensor_py_operators.__rdot__(right, left)
    649 def __rdot__(right, left):
--> 650     return at.math.dense_dot(left, right)

File ~/.conda/envs/dwell-times/lib/python3.10/site-packages/pytensor/tensor/math.py:2106, in dense_dot(a, b)
   2101 a, b = as_tensor_variable(a), as_tensor_variable(b)
   2103 if not isinstance(a.type, DenseTensorType) or not isinstance(
   2104     b.type, DenseTensorType
   2105 ):
-> 2106     raise TypeError("The dense dot product is only supported for dense types")
   2108 if a.ndim == 0 or b.ndim == 0:
   2109     return a * b

TypeError: The dense dot product is only supported for dense types```

ricardoV94 commented 1 year ago

Seems like a bug? I guess gz is sparse in that case and the gradient shouldn't be calling dot? Or it should atleast first convert gz to a dense Tensor.

pymc-devs / pytensor