1""" 2Classes for handling sparse matrices. 3 4To read about different sparse formats, see 5http://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps 6 7""" 8from __future__ import absolute_import, print_function, division 9 10# TODO 11# Automatic methods for determining best sparse format? 12 13import sys 14 15import numpy as np 16from numpy.lib.stride_tricks import as_strided 17from six import integer_types 18from six.moves import xrange 19import scipy.sparse 20 21import theano 22from theano import gof, tensor, scalar, config 23from theano.gradient import DisconnectedType 24from theano.sparse.utils import hash_from_sparse 25from theano.gradient import grad_not_implemented, grad_undefined 26from theano.sparse.type import SparseType, _is_sparse 27 28sparse_formats = ['csc', 'csr'] 29 30 31""" 32Types of sparse matrices to use for testing. 33 34""" 35_mtypes = [scipy.sparse.csc_matrix, scipy.sparse.csr_matrix] 36# _mtypes = [sparse.csc_matrix, sparse.csr_matrix, sparse.dok_matrix, 37# sparse.lil_matrix, sparse.coo_matrix] 38# * new class ``dia_matrix`` : the sparse DIAgonal format 39# * new class ``bsr_matrix`` : the Block CSR format 40_mtype_to_str = {scipy.sparse.csc_matrix: "csc", 41 scipy.sparse.csr_matrix: "csr"} 42 43 44def _is_sparse_variable(x): 45 """ 46 47 Returns 48 ------- 49 boolean 50 True iff x is a L{SparseVariable} (and not a L{tensor.TensorType}, 51 for instance). 52 53 """ 54 if not isinstance(x, gof.Variable): 55 raise NotImplementedError("this function should only be called on " 56 "*variables* (of type sparse.SparseType " 57 "or tensor.TensorType, for instance), not ", 58 x) 59 return isinstance(x.type, SparseType) 60 61 62def _is_dense_variable(x): 63 """ 64 65 Returns 66 ------- 67 boolean 68 True if x is a L{tensor.TensorType} (and not a L{SparseVariable}, 69 for instance). 70 71 """ 72 if not isinstance(x, gof.Variable): 73 raise NotImplementedError("this function should only be called on " 74 "*variables* (of type sparse.SparseType or " 75 "tensor.TensorType, for instance), not ", x) 76 return isinstance(x.type, tensor.TensorType) 77 78 79def _is_dense(x): 80 """ 81 82 Returns 83 ------- 84 boolean 85 True unless x is a L{scipy.sparse.spmatrix} (and not a 86 L{numpy.ndarray}). 87 88 """ 89 if not isinstance(x, (scipy.sparse.spmatrix, np.ndarray)): 90 raise NotImplementedError("this function should only be called on " 91 "sparse.scipy.sparse.spmatrix or " 92 "numpy.ndarray, not,", x) 93 return isinstance(x, np.ndarray) 94 95 96# Wrapper type 97def as_sparse_variable(x, name=None): 98 """ 99 Wrapper around SparseVariable constructor to construct 100 a Variable with a sparse matrix with the same dtype and 101 format. 102 103 Parameters 104 ---------- 105 x 106 A sparse matrix. 107 108 Returns 109 ------- 110 object 111 SparseVariable version of `x`. 112 113 """ 114 115 # TODO 116 # Verify that sp is sufficiently sparse, and raise a 117 # warning if it is not 118 119 if isinstance(x, gof.Apply): 120 if len(x.outputs) != 1: 121 raise ValueError("It is ambiguous which output of a " 122 "multi-output Op has to be fetched.", x) 123 else: 124 x = x.outputs[0] 125 if isinstance(x, gof.Variable): 126 if not isinstance(x.type, SparseType): 127 raise TypeError("Variable type field must be a SparseType.", x, 128 x.type) 129 return x 130 try: 131 return constant(x, name=name) 132 except TypeError: 133 raise TypeError("Cannot convert %s to SparseType" % x, type(x)) 134as_sparse = as_sparse_variable 135 136 137def as_sparse_or_tensor_variable(x, name=None): 138 """ 139 Same as `as_sparse_variable` but if we can't make a 140 sparse variable, we try to make a tensor variable. 141 142 Parameters 143 ---------- 144 x 145 A sparse matrix. 146 147 Returns 148 ------- 149 SparseVariable or TensorVariable version of `x` 150 151 """ 152 153 try: 154 return as_sparse_variable(x, name) 155 except (ValueError, TypeError): 156 return theano.tensor.as_tensor_variable(x, name) 157 158 159def constant(x, name=None): 160 if not isinstance(x, scipy.sparse.spmatrix): 161 raise TypeError("sparse.constant must be called on a " 162 "scipy.sparse.spmatrix") 163 try: 164 return SparseConstant(SparseType(format=x.format, 165 dtype=x.dtype), x.copy(), name=name) 166 except TypeError: 167 raise TypeError("Could not convert %s to SparseType" % x, type(x)) 168 169 170def sp_ones_like(x): 171 """ 172 Construct a sparse matrix of ones with the same sparsity pattern. 173 174 Parameters 175 ---------- 176 x 177 Sparse matrix to take the sparsity pattern. 178 179 Returns 180 ------- 181 A sparse matrix 182 The same as `x` with data changed for ones. 183 184 """ 185 # TODO: don't restrict to CSM formats 186 data, indices, indptr, shape = csm_properties(x) 187 return CSM(format=x.format)(tensor.ones_like(data), indices, indptr, shape) 188 189 190def sp_zeros_like(x): 191 """ 192 Construct a sparse matrix of zeros. 193 194 Parameters 195 ---------- 196 x 197 Sparse matrix to take the shape. 198 199 Returns 200 ------- 201 A sparse matrix 202 The same as `x` with zero entries for all element. 203 204 """ 205 206 # TODO: don't restrict to CSM formats 207 _, _, indptr, shape = csm_properties(x) 208 return CSM(format=x.format)(data=np.array([], dtype=x.type.dtype), 209 indices=np.array([], dtype='int32'), 210 indptr=tensor.zeros_like(indptr), 211 shape=shape) 212 213 214class _sparse_py_operators: 215 T = property(lambda self: transpose(self), 216 doc="Return aliased transpose of self (read-only)") 217 218 def astype(self, dtype): 219 return cast(self, dtype) 220 221 def __neg__(self): 222 return neg(self) 223 224 def __add__(left, right): 225 return add(left, right) 226 227 def __radd__(right, left): 228 return add(left, right) 229 230 def __sub__(left, right): 231 return sub(left, right) 232 233 def __rsub__(right, left): 234 return sub(left, right) 235 236 def __mul__(left, right): 237 return mul(left, right) 238 239 def __rmul__(left, right): 240 return mul(left, right) 241 242 # comparison operators 243 244 def __lt__(self, other): 245 return lt(self, other) 246 247 def __le__(self, other): 248 return le(self, other) 249 250 def __gt__(self, other): 251 return gt(self, other) 252 253 def __ge__(self, other): 254 return ge(self, other) 255 256 # extra pseudo-operator symbols 257 258 def __dot__(left, right): 259 return structured_dot(left, right) 260 261 def __rdot__(right, left): 262 return structured_dot(left, right) 263 264 # N.B. THIS IS COMMENTED OUT ON PURPOSE!!! 265 # Discussion with Fred & James (at least, and maybe others before) 266 # we decided that casting from a sparse to dense should be explicit 267 # because it's usually something you just want to be pretty careful 268 # about, and not to do by accident. 269 # def _as_TensorVariable(self): 270 # return dense_from_sparse(self) 271 272 def toarray(self): 273 return dense_from_sparse(self) 274 shape = property(lambda self: tensor.shape(dense_from_sparse(self))) 275 # don't worry! 276 # the plan is that the ShapeFeature in tensor.opt will do shape propagation 277 # and remove the dense_from_sparse from the graph. This will *NOT* 278 # actually expand your sparse matrix just to get the shape. 279 ndim = property(lambda self: self.type.ndim) 280 dtype = property(lambda self: self.type.dtype) 281 282 # Note that the `size` attribute of sparse matrices behaves differently 283 # from dense matrices: it is the number of elements stored in the matrix 284 # rather than the total number of elements that may be stored. Note also 285 # that stored zeros *do* count in the size. 286 size = property(lambda self: csm_data(self).size) 287 288 def zeros_like(model): 289 return sp_zeros_like(model) 290 291 def __getitem__(self, args): 292 if not isinstance(args, tuple): 293 args = args, 294 295 if len(args) == 2: 296 scalar_arg_1 = (np.isscalar(args[0]) or 297 getattr(args[0], 'type', None) == tensor.iscalar) 298 scalar_arg_2 = (np.isscalar(args[1]) or 299 getattr(args[1], 'type', None) == tensor.iscalar) 300 if scalar_arg_1 and scalar_arg_2: 301 ret = get_item_scalar(self, args) 302 elif isinstance(args[0], list): 303 ret = get_item_2lists(self, args[0], args[1]) 304 else: 305 ret = get_item_2d(self, args) 306 elif isinstance(args[0], list): 307 ret = get_item_list(self, args[0]) 308 else: 309 ret = get_item_2d(self, args) 310 return ret 311 312 313class SparseVariable(_sparse_py_operators, gof.Variable): 314 dtype = property(lambda self: self.type.dtype) 315 format = property(lambda self: self.type.format) 316 317 def __str__(self): 318 return '%s{%s,%s}' % ( 319 self.__class__.__name__, 320 self.format, 321 self.dtype) 322 323 def __repr__(self): 324 return str(self) 325 326 327class SparseConstantSignature(tuple): 328 def __eq__(self, other): 329 (a, b), (x, y) = self, other 330 return (a == x and 331 (b.dtype == y.dtype) and 332 (type(b) == type(y)) and 333 (b.shape == y.shape) and 334 (abs(b - y).sum() < 1e-6 * b.nnz)) 335 336 def __ne__(self, other): 337 return not self == other 338 339 def __hash__(self): 340 (a, b) = self 341 return hash(type(self)) ^ hash(a) ^ hash(type(b)) 342 343 def theano_hash(self): 344 (_, d) = self 345 return hash_from_sparse(d) 346 347 348class SparseConstant(gof.Constant, _sparse_py_operators): 349 dtype = property(lambda self: self.type.dtype) 350 format = property(lambda self: self.type.format) 351 352 def signature(self): 353 assert self.data is not None 354 return SparseConstantSignature((self.type, self.data)) 355 356 def __str__(self): 357 return '%s{%s,%s,shape=%s,nnz=%s}' % ( 358 self.__class__.__name__, 359 self.format, 360 self.dtype, 361 self.data.shape, 362 self.data.nnz) 363 364 def __repr__(self): 365 return str(self) 366 367 368SparseType.Variable = SparseVariable 369SparseType.Constant = SparseConstant 370 371 372# for more dtypes, call SparseType(format, dtype) 373def matrix(format, name=None, dtype=None): 374 if dtype is None: 375 dtype = config.floatX 376 type = SparseType(format=format, dtype=dtype) 377 return type(name) 378 379 380def csc_matrix(name=None, dtype=None): 381 return matrix('csc', name, dtype) 382 383 384def csr_matrix(name=None, dtype=None): 385 return matrix('csr', name, dtype) 386 387 388def bsr_matrix(name=None, dtype=None): 389 return matrix('bsr', name, dtype) 390 391 392# for more dtypes, call SparseType(format, dtype) 393csc_dmatrix = SparseType(format='csc', dtype='float64') 394csr_dmatrix = SparseType(format='csr', dtype='float64') 395bsr_dmatrix = SparseType(format='bsr', dtype='float64') 396csc_fmatrix = SparseType(format='csc', dtype='float32') 397csr_fmatrix = SparseType(format='csr', dtype='float32') 398bsr_fmatrix = SparseType(format='bsr', dtype='float32') 399 400all_dtypes = SparseType.dtype_set 401complex_dtypes = [t for t in all_dtypes if t[:7] == 'complex'] 402float_dtypes = [t for t in all_dtypes if t[:5] == 'float'] 403int_dtypes = [t for t in all_dtypes if t[:3] == 'int'] 404uint_dtypes = [t for t in all_dtypes if t[:4] == 'uint'] 405integer_dtypes = int_dtypes + uint_dtypes 406 407continuous_dtypes = complex_dtypes + float_dtypes 408discrete_dtypes = int_dtypes + uint_dtypes 409 410 411# CONSTRUCTION 412class CSMProperties(gof.Op): 413 # See doc in instance of this Op or function after this class definition. 414 # NOTE 415 # We won't implement infer_shape for this op now. This will 416 # ask that we implement an GetNNZ op, and this op will keep 417 # the dependence on the input of this op. So this won't help 418 # to remove computations in the graph. To remove computation, 419 # we will need to make an infer_sparse_pattern feature to 420 # remove computations. Doing this is trickier then the 421 # infer_shape feature. For example, how do we handle the case 422 # when some op create some 0 values? So there is dependence 423 # on the values themselves. We could write an infer_shape for 424 # the last output that is the shape, but I dough this will 425 # get used. 426 427 # we don't return a view of the shape, we create a new ndarray from the 428 # shape tuple. 429 __props__ = () 430 view_map = {0: [0], 1: [0], 2: [0]} 431 432 """ 433 Indexing to speficied what part of the data parameter 434 should be use to construct the sparse matrix. 435 436 """ 437 438 def __init__(self, kmap=None): 439 if kmap is not None: 440 raise Exception("Do not use kmap, it is removed") 441 442 def make_node(self, csm): 443 csm = as_sparse_variable(csm) 444 assert csm.format in ["csr", "csc"] 445 data = tensor.TensorType(dtype=csm.type.dtype, 446 broadcastable=(False,))() 447 return gof.Apply(self, [csm], 448 [data, tensor.ivector(), 449 tensor.ivector(), tensor.ivector()]) 450 451 def perform(self, node, inputs, out): 452 (csm,) = inputs 453 out[0][0] = csm.data 454 if str(csm.data.dtype) == 'int32': 455 out[0][0] = theano._asarray(out[0][0], dtype='int32') 456 # backport 457 out[1][0] = theano._asarray(csm.indices, dtype='int32') 458 out[2][0] = theano._asarray(csm.indptr, dtype='int32') 459 out[3][0] = theano._asarray(csm.shape, dtype='int32') 460 461 def grad(self, inputs, g): 462 463 # g[1:] is all integers, so their Jacobian in this op 464 # is 0. We thus don't need to worry about what their values 465 # are. 466 467 # if g[0] is disconnected, then this op doesn't contribute 468 # any gradient anywhere. but we know that at least one of 469 # g[1:] is connected, or this grad method wouldn't have been 470 # called, so we should report zeros 471 (csm,) = inputs 472 if isinstance(g[0].type, DisconnectedType): 473 return [csm.zeros_like()] 474 475 data, indices, indptr, shape = csm_properties(csm) 476 return [CSM(csm.format)(g[0], indices, indptr, shape)] 477 478# don't make this a function or it breaks some optimizations below 479csm_properties = CSMProperties() 480""" 481Extract all of .data, .indices, .indptr and .shape field. 482 483For specific field, `csm_data`, `csm_indices`, `csm_indptr` 484and `csm_shape` are provided. 485 486Parameters 487---------- 488csm 489 Sparse matrix in CSR or CSC format. 490 491Returns 492 (data, indices, indptr, shape), the properties of `csm`. 493 494Notes 495----- 496The grad implemented is regular, i.e. not structured. 497`infer_shape` method is not available for this op. 498 499""" 500 501 502def csm_data(csm): 503 """ 504 Return the data field of the sparse variable. 505 506 """ 507 return csm_properties(csm)[0] 508 509 510def csm_indices(csm): 511 """ 512 Return the indices field of the sparse variable. 513 514 """ 515 return csm_properties(csm)[1] 516 517 518def csm_indptr(csm): 519 """ 520 Return the indptr field of the sparse variable. 521 522 """ 523 return csm_properties(csm)[2] 524 525 526def csm_shape(csm): 527 """ 528 Return the shape field of the sparse variable. 529 530 """ 531 return csm_properties(csm)[3] 532 533 534class CSM(gof.Op): 535 # See doc in instance of this Op or function after this class definition. 536 """ 537 Indexing to speficied what part of the data parameter 538 should be used to construct the sparse matrix. 539 540 """ 541 __props__ = ('format',) 542 """ 543 Pre-computed hash value, defined by __init__. 544 545 """ 546 547 def __init__(self, format, kmap=None): 548 if format not in ('csr', 'csc'): 549 raise ValueError("format must be one of: 'csr', 'csc'", format) 550 self.format = format 551 if kmap is not None: 552 raise Exception("Do not use kmap, it is removed") 553 # should view the other inputs too, but viewing multiple 554 # inputs is not currently supported by the destroyhandler 555 self.view_map = {0: [0]} 556 557 def make_node(self, data, indices, indptr, shape): 558 data = tensor.as_tensor_variable(data) 559 560 if not isinstance(indices, gof.Variable): 561 indices_ = np.asarray(indices) 562 indices_32 = theano._asarray(indices, dtype='int32') 563 assert (indices_ == indices_32).all() 564 indices = indices_32 565 if not isinstance(indptr, gof.Variable): 566 indptr_ = np.asarray(indptr) 567 indptr_32 = theano._asarray(indptr, dtype='int32') 568 assert (indptr_ == indptr_32).all() 569 indptr = indptr_32 570 if not isinstance(shape, gof.Variable): 571 shape_ = np.asarray(shape) 572 shape_32 = theano._asarray(shape, dtype='int32') 573 assert (shape_ == shape_32).all() 574 shape = shape_32 575 576 indices = tensor.as_tensor_variable(indices) 577 indptr = tensor.as_tensor_variable(indptr) 578 shape = tensor.as_tensor_variable(shape) 579 580 if data.type.ndim != 1: 581 raise TypeError('data argument must be a vector', data.type, 582 data.type.ndim) 583 if indices.type.ndim != 1 or indices.type.dtype not in discrete_dtypes: 584 raise TypeError('indices must be vector of integers', indices, 585 indices.type) 586 if indptr.type.ndim != 1 or indptr.type.dtype not in discrete_dtypes: 587 raise TypeError('indices must be vector of integers', indptr, 588 indptr.type) 589 if shape.type.ndim != 1 or shape.type.dtype not in discrete_dtypes: 590 raise TypeError('n_rows must be integer type', shape, shape.type) 591 592 return gof.Apply(self, 593 [data, indices, indptr, shape], 594 [SparseType(dtype=data.type.dtype, 595 format=self.format)()]) 596 597 def perform(self, node, inputs, outputs): 598 # for efficiency, if remap does nothing, then do not apply it 599 (data, indices, indptr, shape) = inputs 600 (out,) = outputs 601 602 if len(shape) != 2: 603 raise ValueError('Shape should be an array of length 2') 604 if data.shape != indices.shape: 605 errmsg = ('Data (shape ' + repr(data.shape) + 606 ' must have the same number of elements ' + 607 'as indices (shape' + repr(indices.shape) + 608 ')') 609 raise ValueError(errmsg) 610 if self.format == 'csc': 611 out[0] = scipy.sparse.csc_matrix((data, indices.copy(), 612 indptr.copy()), 613 np.asarray(shape), copy=False) 614 else: 615 assert self.format == 'csr' 616 out[0] = scipy.sparse.csr_matrix((data, indices.copy(), 617 indptr.copy()), shape.copy(), 618 copy=False) 619 620 def connection_pattern(self, node): 621 return [[True], [False], [False], [False]] 622 623 def grad(self, inputs, gout): 624 (x_data, x_indices, x_indptr, x_shape) = inputs 625 (g_out,) = gout 626 g_data, g_indices, g_indptr, g_shape = csm_properties(g_out) 627 # unpack the data vector and wrap it as a 1d TensorType 628 g_data = csm_grad()(x_data, x_indices, x_indptr, x_shape, 629 g_data, g_indices, g_indptr, g_shape) 630 return [g_data, DisconnectedType()(), DisconnectedType()(), DisconnectedType()()] 631 632 def infer_shape(self, node, shapes): 633 # node.inputs[3] is of length as we only support sparse matrix. 634 return [(node.inputs[3][0], node.inputs[3][1])] 635 636CSC = CSM('csc') 637""" 638Construct a CSC matrix from the internal representation. 639 640Parameters 641---------- 642data 643 One dimensional tensor representing the data of the sparse matrix to 644 construct. 645indices 646 One dimensional tensor of integers representing the indices of the sparse 647 matrix to construct. 648indptr 649 One dimensional tensor of integers representing the indice pointer for 650 the sparse matrix to construct. 651shape 652 One dimensional tensor of integers representing the shape of the sparse 653 matrix to construct. 654 655Returns 656------- 657sparse matrix 658 A sparse matrix having the properties specified by the inputs. 659 660Notes 661----- 662The grad method returns a dense vector, so it provides a regular grad. 663 664""" 665 666CSR = CSM('csr') 667""" 668Construct a CSR matrix from the internal representation. 669 670Parameters 671---------- 672data 673 One dimensional tensor representing the data of the sparse matrix to 674 construct. 675indices 676 One dimensional tensor of integers representing the indices of the sparse 677 matrix to construct. 678indptr 679 One dimensional tensor of integers representing the indice pointer for 680 the sparse matrix to construct. 681shape 682 One dimensional tensor of integers representing the shape of the sparse 683 matrix to construct. 684 685Returns 686------- 687sparse matrix 688 A sparse matrix having the properties specified by the inputs. 689 690Notes 691----- 692The grad method returns a dense vector, so it provides a regular grad. 693 694""" 695 696 697class CSMGrad(gof.op.Op): 698 # Note 699 # This Op computes the gradient of the CSM Op. CSM creates a matrix from 700 # data, indices, and indptr vectors; it's gradient is the gradient of 701 # the data vector only. There are two complexities to calculate this 702 # gradient: 703 # 1. The gradient may be sparser than the input matrix defined by (data, 704 # indices, indptr). In this case, the data vector of the gradient will have 705 # less elements than the data vector of the input because sparse formats 706 # remove 0s. Since we are only returning the gradient of the data vector, 707 # the relevant 0s need to be added back. 708 # 2. The elements in the sparse dimension are not guaranteed to be sorted. 709 # Therefore, the input data vector may have a different order than the 710 # gradient data vector. 711 __props__ = () 712 713 def __init__(self, kmap=None): 714 if kmap is not None: 715 raise Exception("Do not use kmap, it is removed") 716 # This class always allocate a new output. 717 # I keep this here to help GD understand what this kmap think is. 718 # if self.kmap is None: 719 # self.view_map = {0: [1]} 720 721 def make_node(self, x_data, x_indices, x_indptr, x_shape, 722 g_data, g_indices, g_indptr, g_shape): 723 gout_data = g_data.type() 724 return gof.Apply(self, [x_data, x_indices, x_indptr, x_shape, 725 g_data, g_indices, g_indptr, g_shape], [gout_data]) 726 727 def perform(self, node, inputs, outputs): 728 (x_data, x_indices, x_indptr, x_shape, 729 g_data, g_indices, g_indptr, g_shape) = inputs 730 (g_out,) = outputs 731 if len(x_indptr) - 1 == x_shape[0]: 732 sp_dim = x_shape[1] 733 else: 734 sp_dim = x_shape[0] 735 736 g_row = np.zeros(sp_dim, dtype=g_data.dtype) 737 gout_data = np.zeros(x_data.shape, dtype=node.outputs[0].dtype) 738 739 for i in range(len(x_indptr) - 1): 740 for j_ptr in range(g_indptr[i], g_indptr[i + 1]): 741 g_row[g_indices[j_ptr]] += g_data[j_ptr] 742 743 for j_ptr in range(x_indptr[i], x_indptr[i + 1]): 744 gout_data[j_ptr] = g_row[x_indices[j_ptr]] 745 746 for j_ptr in range(g_indptr[i], g_indptr[i + 1]): 747 g_row[g_indices[j_ptr]] = 0 748 749 g_out[0] = gout_data 750 751 def infer_shape(self, node, shapes): 752 return [shapes[1]] 753 754csm_grad = CSMGrad 755 756 757class Cast(gof.op.Op): 758 # See doc in instance of this Op or function after this class definition. 759 __props__ = ("out_type",) 760 761 def __init__(self, out_type): 762 self.out_type = out_type 763 764 def make_node(self, x): 765 x = as_sparse_variable(x) 766 assert x.format in ["csr", "csc"] 767 return gof.Apply( 768 self, [x], 769 [SparseType(dtype=self.out_type, format=x.format)()]) 770 771 def perform(self, node, inputs, outputs): 772 (x,) = inputs 773 (out,) = outputs 774 assert _is_sparse(x) 775 out[0] = x.astype(self.out_type) 776 777 def grad(self, inputs, outputs_gradients): 778 gz = outputs_gradients[0] 779 780 if gz.dtype in complex_dtypes: 781 raise NotImplementedError("grad not implemented for complex types") 782 if inputs[0].dtype in complex_dtypes: 783 raise NotImplementedError("grad not implemented for complex types") 784 785 if gz.dtype in discrete_dtypes: 786 if inputs[0].dtype in discrete_dtypes: 787 return [inputs[0].zeros_like(dtype=theano.config.floatX)] 788 else: 789 return [inputs[0].zeros_like()] 790 else: 791 if inputs[0].dtype in discrete_dtypes: 792 return [gz] 793 else: 794 return [Cast(inputs[0].dtype)(gz)] 795 796 def infer_shape(self, node, ins_shapes): 797 return ins_shapes 798 799 def __str__(self): 800 return "%s(%s)" % (self.__class__.__name__, self.out_type) 801 802bcast = Cast('int8') 803wcast = Cast('int16') 804icast = Cast('int32') 805lcast = Cast('int64') 806fcast = Cast('float32') 807dcast = Cast('float64') 808ccast = Cast('complex64') 809zcast = Cast('complex128') 810 811 812def cast(variable, dtype): 813 """ 814 Cast sparse variable to the desired dtype. 815 816 Parameters 817 ---------- 818 variable 819 Sparse matrix. 820 dtype 821 The dtype wanted. 822 823 Returns 824 ------- 825 Same as `x` but having `dtype` as dtype. 826 827 Notes 828 ----- 829 The grad implemented is regular, i.e. not structured. 830 831 """ 832 return Cast(dtype)(variable) 833 834# 835# Conversion 836# 837 838 839class DenseFromSparse(gof.op.Op): 840 # See doc in instance of this Op or function after this class definition. 841 __props__ = () # We don't put sparse_grad in the props. 842 843 def __init__(self, structured=True): 844 self.sparse_grad = structured 845 846 def __str__(self): 847 return "%s{structured_grad=%s}" % ( 848 self.__class__.__name__, 849 self.sparse_grad) 850 851 def make_node(self, x): 852 x = as_sparse_variable(x) 853 return gof.Apply(self, 854 [x], 855 [tensor.TensorType(dtype=x.type.dtype, 856 broadcastable=(False, False))()]) 857 858 def perform(self, node, inputs, outputs): 859 (x,) = inputs 860 (out,) = outputs 861 if _is_dense(x): 862 print(( 863 "WARNING: You just called DenseFromSparse on a dense matrix." 864 ), file=sys.stderr) 865 out[0] = x 866 else: 867 out[0] = x.toarray() 868 assert _is_dense(out[0]) 869 870 def grad(self, inputs, gout): 871 (x,) = inputs 872 (gz,) = gout 873 if self.sparse_grad: 874 left = sp_ones_like(x) 875 right = gz 876 877 # Do upcasting if necessary to avoid an unimplemented case 878 # of mul 879 880 if right.dtype == 'float64' and left.dtype == 'float32': 881 left = left.astype('float64') 882 883 if right.dtype == 'float32' and left.dtype == 'float64': 884 right = right.astype('float64') 885 886 return [left * right] 887 else: 888 return [SparseFromDense(x.type.format)(gz)] 889 890 def infer_shape(self, node, shapes): 891 return [shapes[0]] 892 893dense_from_sparse = DenseFromSparse() 894""" 895Convert a sparse matrix to a dense one. 896 897Parameters 898---------- 899x 900 A sparse matrix. 901 902Returns 903------- 904theano.tensor.matrix 905 A dense matrix, the same as `x`. 906 907Notes 908----- 909The grad implementation can be controlled through the constructor via the 910`structured` parameter. `True` will provide a structured grad while `False` 911will provide a regular grad. By default, the grad is structured. 912 913""" 914 915 916class SparseFromDense(gof.op.Op): 917 918 __props__ = () 919 920 def __init__(self, format): 921 self.format = format 922 923 def __str__(self): 924 return "%s{%s}" % ( 925 self.__class__.__name__, 926 self.format) 927 928 def make_node(self, x): 929 x = tensor.as_tensor_variable(x) 930 if x.ndim > 2: 931 raise TypeError( 932 "Theano does not have sparse tensor types with more " 933 "than 2 dimensions, but %s.ndim = %i" % (x, x.ndim)) 934 elif x.ndim == 1: 935 x = x.dimshuffle('x', 0) 936 elif x.ndim == 0: 937 x = x.dimshuffle('x', 'x') 938 else: 939 assert x.ndim == 2 940 941 return gof.Apply(self, 942 [x], 943 [SparseType(dtype=x.type.dtype, 944 format=self.format)()]) 945 946 def perform(self, node, inputs, outputs): 947 (x,) = inputs 948 (out,) = outputs 949 out[0] = SparseType.format_cls[self.format](x) 950 951 def grad(self, inputs, gout): 952 (x,) = inputs 953 (gz,) = gout 954 gx = dense_from_sparse(gz) 955 gx = tensor.patternbroadcast(gx, x.broadcastable) 956 return gx, 957 958 def infer_shape(self, node, shapes): 959 return [shapes[0]] 960 961csr_from_dense = SparseFromDense('csr') 962""" 963Convert a dense matrix to a sparse csr matrix. 964 965Parameters 966---------- 967x 968 A dense matrix. 969 970Returns 971------- 972sparse matrix 973 The same as `x` in a sparse csr matrix format. 974 975""" 976 977csc_from_dense = SparseFromDense('csc') 978""" 979Convert a dense matrix to a sparse csc matrix. 980 981Parameters 982---------- 983x 984 A dense matrix. 985 986Returns 987------- 988sparse matrix 989 The same as `x` in a sparse csc matrix format. 990 991""" 992 993 994# Indexing 995class GetItemList(gof.op.Op): 996 997 __props__ = () 998 999 def infer_shape(self, node, shapes): 1000 return [(shapes[1][0], shapes[0][1])] 1001 1002 def make_node(self, x, index): 1003 x = as_sparse_variable(x) 1004 assert x.format in ["csr", "csc"] 1005 1006 ind = tensor.as_tensor_variable(index) 1007 assert ind.ndim == 1 1008 assert ind.dtype in integer_dtypes 1009 1010 return gof.Apply(self, [x, ind], [x.type()]) 1011 1012 def perform(self, node, inp, outputs): 1013 (out,) = outputs 1014 x = inp[0] 1015 indices = inp[1] 1016 assert _is_sparse(x) 1017 out[0] = x[indices] 1018 1019 def grad(self, inputs, g_outputs): 1020 x, indices = inputs 1021 gout, = g_outputs 1022 return [get_item_list_grad(x, indices, gout), 1023 grad_undefined(self, 1, indices, "No gradient for this input")] 1024 1025get_item_list = GetItemList() 1026""" 1027Select row of sparse matrix, returning them as a new sparse matrix. 1028 1029Parameters 1030---------- 1031x 1032 Sparse matrix. 1033index 1034 List of rows. 1035 1036Returns 1037------- 1038sparse matrix 1039 The corresponding rows in `x`. 1040 1041""" 1042 1043 1044class GetItemListGrad(gof.op.Op): 1045 1046 __props__ = () 1047 1048 def infer_shape(self, node, shapes): 1049 return [(shapes[0])] 1050 1051 def make_node(self, x, index, gz): 1052 x = as_sparse_variable(x) 1053 gz = as_sparse_variable(gz) 1054 1055 assert x.format in ["csr", "csc"] 1056 assert gz.format in ["csr", "csc"] 1057 1058 ind = tensor.as_tensor_variable(index) 1059 assert ind.ndim == 1 1060 assert ind.dtype in integer_dtypes 1061 1062 scipy_ver = [int(n) for n in scipy.__version__.split('.')[:2]] 1063 1064 if not scipy_ver >= [0, 13]: 1065 raise NotImplementedError("Scipy version is to old") 1066 1067 return gof.Apply(self, [x, ind, gz], [x.type()]) 1068 1069 def perform(self, node, inp, outputs): 1070 (out,) = outputs 1071 x = inp[0] 1072 indices = inp[1] 1073 gz = inp[2] 1074 1075 if x.format in ["csr"]: 1076 y = scipy.sparse.csr_matrix((x.shape[0], x.shape[1])) 1077 else: 1078 y = scipy.sparse.csc_matrix((x.shape[0], x.shape[1])) 1079 for a in range(0, len(indices)): 1080 y[indices[a]] = gz[a] 1081 1082 out[0] = y 1083 1084get_item_list_grad = GetItemListGrad() 1085 1086 1087class GetItem2Lists(gof.op.Op): 1088 1089 __props__ = () 1090 1091 def make_node(self, x, ind1, ind2): 1092 x = as_sparse_variable(x) 1093 assert x.format in ["csr", "csc"] 1094 ind1 = tensor.as_tensor_variable(ind1) 1095 ind2 = tensor.as_tensor_variable(ind2) 1096 assert ind1.dtype in integer_dtypes 1097 assert ind2.dtype in integer_dtypes 1098 1099 return gof.Apply(self, [x, ind1, ind2], 1100 [theano.tensor.vector()]) 1101 1102 def perform(self, node, inp, outputs): 1103 (out,) = outputs 1104 x = inp[0] 1105 ind1 = inp[1] 1106 ind2 = inp[2] 1107 out[0] = np.asarray(x[ind1, ind2]).flatten() 1108 """ 1109 Here scipy returns the corresponding elements in a matrix which isn't 1110 what we are aiming for. Using asarray and flatten, out[0] becomes an 1111 array. 1112 1113 """ 1114 def grad(self, inputs, g_outputs): 1115 x, ind1, ind2 = inputs 1116 gout, = g_outputs 1117 return [get_item_2lists_grad(x, ind1, ind2, gout), 1118 grad_undefined(self, 1, ind1, "No gradient for this input"), 1119 grad_undefined(self, 1, ind2, "No gradient for this input")] 1120 1121get_item_2lists = GetItem2Lists() 1122""" 1123Select elements of sparse matrix, returning them in a vector. 1124 1125Parameters 1126---------- 1127x 1128 Sparse matrix. 1129index 1130 List of two lists, first list indicating the row of each element and second 1131 list indicating its column. 1132 1133Returns 1134------- 1135theano.tensor.vector 1136 The corresponding elements in `x`. 1137 1138""" 1139 1140 1141class GetItem2ListsGrad(gof.op.Op): 1142 1143 __props__ = () 1144 1145 def infer_shape(self, node, shapes): 1146 return [(shapes[0])] 1147 1148 def make_node(self, x, ind1, ind2, gz): 1149 x = as_sparse_variable(x) 1150 1151 assert x.format in ["csr", "csc"] 1152 1153 ind1 = tensor.as_tensor_variable(ind1) 1154 ind2 = tensor.as_tensor_variable(ind2) 1155 assert ind1.ndim == 1 1156 assert ind2.ndim == 1 1157 assert ind1.dtype in integer_dtypes 1158 assert ind2.dtype in integer_dtypes 1159 1160 return gof.Apply(self, [x, ind1, ind2, gz], [x.type()]) 1161 1162 def perform(self, node, inp, outputs): 1163 (out,) = outputs 1164 x = inp[0] 1165 ind1 = inp[1] 1166 ind2 = inp[2] 1167 gz = inp[3] 1168 1169 if x.format in ["csr"]: 1170 y = scipy.sparse.csr_matrix((x.shape[0], x.shape[1])) 1171 else: 1172 y = scipy.sparse.csc_matrix((x.shape[0], x.shape[1])) 1173 z = 0 1174 for z in range(0, len(ind1)): 1175 y[(ind1[z], ind2[z])] = gz[z] 1176 1177 out[0] = y 1178 1179get_item_2lists_grad = GetItem2ListsGrad() 1180 1181 1182class GetItem2d(gof.op.Op): 1183 # See doc in instance of this Op or function after this class definition. 1184 1185 __props__ = () 1186 1187# Fred:Too complicated for now. If you need it, look at 1188# the Subtensor.infer_shape. 1189# def infer_shape(self, node, i0_shapes): 1190# return i0_shapes 1191 def make_node(self, x, index): 1192 scipy_ver = [int(n) for n in scipy.__version__.split('.')[:2]] 1193 x = as_sparse_variable(x) 1194 assert x.format in ["csr", "csc"] 1195 assert len(index) in [1, 2] 1196 1197 input_op = [x] 1198 generic_None = theano.gof.Constant(theano.gof.generic, None) 1199 1200 for ind in index: 1201 if isinstance(ind, slice): 1202 # in case of slice is written in theano variable 1203 start = ind.start 1204 stop = ind.stop 1205 step = ind.step 1206 # If start or stop or step are None, make them a Generic 1207 # constant. Else, they should be converted to Tensor Variables 1208 # of dimension 1 and int/uint dtype. 1209 if scipy_ver < [0, 14] and ind.step is not None: 1210 raise ValueError( 1211 'Slice with step is not support with current' 1212 ' version of Scipy.') 1213 if ind.step is None or ind.step == 1: 1214 step = generic_None 1215 else: 1216 if not isinstance(step, gof.Variable): 1217 step = tensor.as_tensor_variable(step) 1218 if not (step.ndim == 0 and step.dtype in 1219 tensor.discrete_dtypes): 1220 raise ValueError(( 1221 "Impossible to index into a sparse matrix with " 1222 "slice where step=%s" % step), 1223 step.ndim, step.dtype) 1224 1225 if start is None: 1226 start = generic_None 1227 else: 1228 if not isinstance(start, gof.Variable): 1229 start = tensor.as_tensor_variable(start) 1230 if not (start.ndim == 0 and start.dtype in 1231 tensor.discrete_dtypes): 1232 raise ValueError(( 1233 "Impossible to index into a sparse matrix with " 1234 "slice where start=%s" % start), 1235 start.ndim, start.dtype) 1236 1237 if stop is None: 1238 stop = generic_None 1239 else: 1240 if not isinstance(stop, gof.Variable): 1241 stop = tensor.as_tensor_variable(stop) 1242 if not (stop.ndim == 0 and stop.dtype in 1243 tensor.discrete_dtypes): 1244 raise ValueError(( 1245 "Impossible to index into a sparse matrix with " 1246 "slice where stop=%s" % stop), 1247 stop.ndim, stop.dtype) 1248 1249 elif ((isinstance(ind, gof.Variable) and 1250 getattr(ind, 'ndim', -1) == 0) or 1251 np.isscalar(ind)): 1252 raise NotImplementedError( 1253 'Theano has no sparse vector' + 1254 'Use X[a:b, c:d], X[a:b, c:c+1] or X[a:b] instead.') 1255 else: 1256 raise ValueError(( 1257 'Advanced indexing is not implemented for sparse ' 1258 'matrices. Argument not supported: %s' % ind)) 1259 input_op += [start, stop, step] 1260 if len(index) == 1: 1261 input_op += [generic_None, generic_None, generic_None] 1262 1263 return gof.Apply(self, input_op, [x.type()]) 1264 1265 def perform(self, node, inputs, outputs): 1266 (x, start1, stop1, step1, start2, stop2, step2) = inputs 1267 (out,) = outputs 1268 assert _is_sparse(x) 1269 out[0] = x[start1:stop1:step1, start2:stop2:step2] 1270 1271get_item_2d = GetItem2d() 1272""" 1273Implement a subtensor of sparse variable, returning a sparse matrix. 1274 1275If you want to take only one element of a sparse matrix see 1276`GetItemScalar` that returns a tensor scalar. 1277 1278Parameters 1279---------- 1280x 1281 Sparse matrix. 1282index 1283 Tuple of slice object. 1284 1285Returns 1286------- 1287sparse matrix 1288 The corresponding slice in `x`. 1289 1290 1291Notes 1292----- 1293Subtensor selection always returns a matrix, so indexing with [a:b, c:d] 1294is forced. If one index is a scalar, for instance, x[a:b, c] or x[a, b:c], 1295an error will be raised. Use instead x[a:b, c:c+1] or x[a:a+1, b:c]. 1296 1297The above indexing methods are not supported because the return value 1298would be a sparse matrix rather than a sparse vector, which is a 1299deviation from numpy indexing rule. This decision is made largely 1300to preserve consistency between numpy and theano. This may be revised 1301when sparse vectors are supported. 1302 1303The grad is not implemented for this op. 1304 1305""" 1306 1307 1308class GetItemScalar(gof.op.Op): 1309 # See doc in instance of this Op or function after this class definition. 1310 __props__ = () 1311 1312 def infer_shape(self, node, shapes): 1313 return [()] 1314 1315 def make_node(self, x, index): 1316 x = as_sparse_variable(x) 1317 assert x.format in ["csr", "csc"] 1318 assert len(index) == 2 1319 1320 input_op = [x] 1321 1322 for ind in index: 1323 1324 if isinstance(ind, slice): 1325 raise Exception("GetItemScalar called with a slice as index!") 1326 1327 # in case of indexing using int instead of theano variable 1328 elif isinstance(ind, integer_types): 1329 ind = theano.tensor.constant(ind) 1330 input_op += [ind] 1331 1332 # in case of indexing using theano variable 1333 elif ind.ndim == 0: 1334 input_op += [ind] 1335 else: 1336 raise NotImplementedError 1337 1338 return gof.Apply(self, input_op, [tensor.scalar(dtype=x.dtype)]) 1339 1340 def perform(self, node, inputs, outputs): 1341 (x, ind1, ind2) = inputs 1342 (out,) = outputs 1343 assert _is_sparse(x) 1344 out[0] = theano._asarray(x[ind1, ind2], x.dtype) 1345 1346get_item_scalar = GetItemScalar() 1347""" 1348Implement a subtensor of a sparse variable that takes two scalars as index and 1349returns a scalar. 1350 1351If you want to take a slice of a sparse matrix see `GetItem2d` that returns a 1352sparse matrix. 1353 1354Parameters 1355---------- 1356x 1357 Sparse matrix. 1358index 1359 Tuple of scalars. 1360 1361Returns 1362------- 1363theano.tensor.scalar 1364 The corresponding item in `x`. 1365 1366Notes 1367----- 1368The grad is not implemented for this op. 1369 1370""" 1371 1372 1373# Linear Algebra 1374class Transpose(gof.op.Op): 1375 # See doc in instance of this Op or function after this class definition. 1376 view_map = {0: [0]} 1377 1378 format_map = {'csr': 'csc', 1379 'csc': 'csr'} 1380 __props__ = () 1381 1382 def __str__(self): 1383 return "Sparse" + self.__class__.__name__ 1384 1385 def make_node(self, x): 1386 x = as_sparse_variable(x) 1387 assert x.format in ["csr", "csc"] 1388 return gof.Apply(self, 1389 [x], 1390 [SparseType(dtype=x.type.dtype, 1391 format=self.format_map[x.type.format])()]) 1392 1393 def perform(self, node, inputs, outputs): 1394 (x,) = inputs 1395 (out,) = outputs 1396 assert _is_sparse(x) 1397 out[0] = x.transpose() 1398 1399 def grad(self, inputs, gout): 1400 (x,) = inputs 1401 (gz,) = gout 1402 assert _is_sparse_variable(x) and _is_sparse_variable(gz) 1403 return transpose(gz), 1404 1405 def infer_shape(self, node, shapes): 1406 return [shapes[0][::-1]] 1407transpose = Transpose() 1408""" 1409Return the transpose of the sparse matrix. 1410 1411Parameters 1412---------- 1413x 1414 Sparse matrix. 1415 1416Returns 1417------- 1418sparse matrix 1419 `x` transposed. 1420 1421Notes 1422----- 1423The returned matrix will not be in the same format. `csc` matrix will be changed 1424in `csr` matrix and `csr` matrix in `csc` matrix. 1425 1426The grad is regular, i.e. not structured. 1427 1428""" 1429 1430 1431class Neg(gof.op.Op): 1432 # See doc in instance of this Op or function after this class definition. 1433 1434 __props__ = () 1435 1436 def __str__(self): 1437 return "Sparse" + self.__class__.__name__ 1438 1439 def make_node(self, x): 1440 x = as_sparse_variable(x) 1441 assert x.format in ["csr", "csc"] 1442 return gof.Apply(self, [x], [x.type()]) 1443 1444 def perform(self, node, inputs, outputs): 1445 (x,) = inputs 1446 (out,) = outputs 1447 assert _is_sparse(x) 1448 out[0] = -x 1449 1450 def grad(self, inputs, gout): 1451 (x,) = inputs 1452 (gz,) = gout 1453 assert _is_sparse_variable(x) and _is_sparse_variable(gz) 1454 return -gz, 1455 1456 def infer_shape(self, node, shapes): 1457 return [shapes[0]] 1458neg = Neg() 1459""" 1460Return the negation of the sparse matrix. 1461 1462Parameters 1463---------- 1464x 1465 Sparse matrix. 1466 1467Returns 1468------- 1469sparse matrix 1470 -`x`. 1471 1472Notes 1473----- 1474The grad is regular, i.e. not structured. 1475 1476""" 1477 1478 1479class ColScaleCSC(gof.op.Op): 1480 # Scale each columns of a sparse matrix by the corresponding 1481 # element of a dense vector 1482 1483 # :param x: A sparse matrix. 1484 # :param s: A dense vector with length equal to the number 1485 # of columns of `x`. 1486 1487 # :return: A sparse matrix in the same format as `x` which 1488 # each column had been multiply by the corresponding 1489 # element of `s`. 1490 1491 # :note: The grad implemented is structured. 1492 1493 __props__ = () 1494 1495 def make_node(self, x, s): 1496 if x.format != 'csc': 1497 raise ValueError('x was not a csc matrix') 1498 return gof.Apply(self, [x, s], [x.type()]) 1499 1500 def perform(self, node, inputs, outputs): 1501 (x, s) = inputs 1502 (z,) = outputs 1503 M, N = x.shape 1504 assert x.format == 'csc' 1505 assert s.shape == (N, ) 1506 1507 y = x.copy() 1508 1509 for j in xrange(0, N): 1510 y.data[y.indptr[j]: y.indptr[j + 1]] *= s[j] 1511 1512 z[0] = y 1513 1514 def grad(self, inputs, gout): 1515 (x, s) = inputs 1516 (gz,) = gout 1517 return [col_scale(gz, s), sp_sum(x * gz, axis=0)] 1518 1519 def infer_shape(self, node, ins_shapes): 1520 return [ins_shapes[0]] 1521 1522 1523class RowScaleCSC(gof.op.Op): 1524 # Scale each row of a sparse matrix by the corresponding element of 1525 # a dense vector 1526 1527 # :param x: A sparse matrix. 1528 # :param s: A dense vector with length equal to the number 1529 # of rows of `x`. 1530 1531 # :return: A sparse matrix in the same format as `x` which 1532 # each row had been multiply by the corresponding 1533 # element of `s`. 1534 1535 # :note: The grad implemented is structured. 1536 1537 view_map = {0: [0]} 1538 __props__ = () 1539 1540 def make_node(self, x, s): 1541 x = as_sparse_variable(x) 1542 assert x.format in ["csr", "csc"] 1543 return gof.Apply(self, [x, s], [x.type()]) 1544 1545 def perform(self, node, inputs, outputs): 1546 (x, s) = inputs 1547 (z,) = outputs 1548 M, N = x.shape 1549 assert x.format == 'csc' 1550 assert s.shape == (M,) 1551 1552 indices = x.indices 1553 indptr = x.indptr 1554 1555 y_data = x.data.copy() 1556 1557 for j in xrange(0, N): 1558 for i_idx in xrange(indptr[j], indptr[j + 1]): 1559 y_data[i_idx] *= s[indices[i_idx]] 1560 1561 z[0] = scipy.sparse.csc_matrix((y_data, indices, indptr), (M, N)) 1562 1563 def grad(self, inputs, gout): 1564 (x, s) = inputs 1565 (gz,) = gout 1566 return [row_scale(gz, s), sp_sum(x * gz, axis=1)] 1567 1568 def infer_shape(self, node, ins_shapes): 1569 return [ins_shapes[0]] 1570 1571 1572def col_scale(x, s): 1573 """ 1574 Scale each columns of a sparse matrix by the corresponding element of a 1575 dense vector. 1576 1577 Parameters 1578 ---------- 1579 x 1580 A sparse matrix. 1581 s 1582 A dense vector with length equal to the number of columns of `x`. 1583 1584 Returns 1585 ------- 1586 A sparse matrix in the same format as `x` which each column had been 1587 multiply by the corresponding element of `s`. 1588 1589 Notes 1590 ----- 1591 The grad implemented is structured. 1592 1593 """ 1594 1595 if x.format == 'csc': 1596 return ColScaleCSC()(x, s) 1597 elif x.format == 'csr': 1598 return RowScaleCSC()(x.T, s).T 1599 else: 1600 raise NotImplementedError() 1601 1602 1603def row_scale(x, s): 1604 """ 1605 Scale each row of a sparse matrix by the corresponding element of 1606 a dense vector. 1607 1608 Parameters 1609 ---------- 1610 x 1611 A sparse matrix. 1612 s 1613 A dense vector with length equal to the number of rows of `x`. 1614 1615 Returns 1616 ------- 1617 A sparse matrix 1618 A sparse matrix in the same format as `x` whose each row has been 1619 multiplied by the corresponding element of `s`. 1620 1621 Notes 1622 ----- 1623 The grad implemented is structured. 1624 1625 """ 1626 return col_scale(x.T, s).T 1627 1628 1629class SpSum(gof.op.Op): 1630 # See doc in instance of this Op or function after this class definition. 1631 1632 __props__ = ("axis",) 1633 # WARNING: judgement call... 1634 # We are not using the structured in the comparison or hashing 1635 # because it doesn't change the perform method therefore, we 1636 # *do* want Sums with different structured values to be merged 1637 # by the merge optimization and this requires them to compare equal. 1638 1639 def __init__(self, axis=None, sparse_grad=True): 1640 super(SpSum, self).__init__() 1641 self.axis = axis 1642 self.structured = sparse_grad 1643 if self.axis not in (None, 0, 1): 1644 raise ValueError('Illegal value for self.axis.') 1645 1646 def make_node(self, x): 1647 x = as_sparse_variable(x) 1648 assert x.format in ["csr", "csc"] 1649 b = () 1650 if self.axis is not None: 1651 b = (False,) 1652 1653 z = tensor.TensorType(broadcastable=b, dtype=x.dtype)() 1654 return gof.Apply(self, [x], [z]) 1655 1656 def perform(self, node, inputs, outputs): 1657 (x,) = inputs 1658 (z,) = outputs 1659 if self.axis is None: 1660 z[0] = np.asarray(x.sum()) 1661 else: 1662 z[0] = np.asarray(x.sum(self.axis)).ravel() 1663 1664 def grad(self, inputs, gout): 1665 (x,) = inputs 1666 (gz,) = gout 1667 if x.dtype not in continuous_dtypes: 1668 return [x.zeros_like(dtype=theano.config.floatX)] 1669 if self.structured: 1670 if self.axis is None: 1671 r = gz * theano.sparse.sp_ones_like(x) 1672 elif self.axis == 0: 1673 r = col_scale(theano.sparse.sp_ones_like(x), gz) 1674 elif self.axis == 1: 1675 r = row_scale(theano.sparse.sp_ones_like(x), gz) 1676 else: 1677 raise ValueError('Illegal value for self.axis.') 1678 else: 1679 o_format = x.format 1680 x = dense_from_sparse(x) 1681 if _is_sparse_variable(gz): 1682 gz = dense_from_sparse(gz) 1683 if self.axis is None: 1684 r = tensor.second(x, gz) 1685 else: 1686 ones = tensor.ones_like(x) 1687 if self.axis == 0: 1688 r = tensor.addbroadcast(gz.dimshuffle('x', 0), 0) * ones 1689 elif self.axis == 1: 1690 r = tensor.addbroadcast(gz.dimshuffle(0, 'x'), 1) * ones 1691 else: 1692 raise ValueError('Illegal value for self.axis.') 1693 r = SparseFromDense(o_format)(r) 1694 return [r] 1695 1696 def infer_shape(self, node, shapes): 1697 r = None 1698 if self.axis is None: 1699 r = [()] 1700 elif self.axis == 0: 1701 r = [(shapes[0][1],)] 1702 else: 1703 r = [(shapes[0][0],)] 1704 return r 1705 1706 def __str__(self): 1707 return self.__class__.__name__ + "{axis=%s}" % str(self.axis) 1708 1709 1710def sp_sum(x, axis=None, sparse_grad=False): 1711 """ 1712 Calculate the sum of a sparse matrix along the specified axis. 1713 1714 It operates a reduction along the specified axis. When `axis` is `None`, 1715 it is applied along all axes. 1716 1717 Parameters 1718 ---------- 1719 x 1720 Sparse matrix. 1721 axis 1722 Axis along which the sum is applied. Integer or `None`. 1723 sparse_grad : bool 1724 `True` to have a structured grad. 1725 1726 Returns 1727 ------- 1728 object 1729 The sum of `x` in a dense format. 1730 1731 Notes 1732 ----- 1733 The grad implementation is controlled with the `sparse_grad` parameter. 1734 `True` will provide a structured grad and `False` will provide a regular 1735 grad. For both choices, the grad returns a sparse matrix having the same 1736 format as `x`. 1737 1738 This op does not return a sparse matrix, but a dense tensor matrix. 1739 1740 """ 1741 1742 return SpSum(axis, sparse_grad)(x) 1743 1744 1745class Diag(gof.op.Op): 1746 # See doc in instance of this Op or function after this class definition. 1747 __props__ = () 1748 1749 def make_node(self, x): 1750 x = as_sparse_variable(x) 1751 assert x.format in ["csr", "csc"] 1752 return gof.Apply(self, [x], [tensor.tensor(broadcastable=(False,), 1753 dtype=x.dtype)]) 1754 1755 def perform(self, node, inputs, outputs): 1756 (x,) = inputs 1757 (z,) = outputs 1758 N, M = x.shape 1759 if N != M: 1760 raise ValueError('Diag only apply on square matrix') 1761 z[0] = x.diagonal() 1762 1763 def grad(self, inputs, gout): 1764 (x,) = inputs 1765 (gz,) = gout 1766 return [square_diagonal(gz)] 1767 1768 def infer_shape(self, nodes, shapes): 1769 return [(tensor.minimum(*shapes[0]), )] 1770 1771diag = Diag() 1772""" 1773Extract the diagonal of a square sparse matrix as a dense vector. 1774 1775Parameters 1776---------- 1777x 1778 A square sparse matrix in csc format. 1779 1780Returns 1781------- 1782theano.tensor.vector 1783 A dense vector representing the diagonal elements. 1784 1785Notes 1786----- 1787The grad implemented is regular, i.e. not structured, since the output is a 1788dense vector. 1789 1790""" 1791 1792 1793class SquareDiagonal(gof.op.Op): 1794 # See doc in instance of this Op or function after this class definition. 1795 1796 __props__ = () 1797 1798 def make_node(self, diag): 1799 diag = tensor.as_tensor_variable(diag) 1800 if diag.type.ndim != 1: 1801 raise TypeError('data argument must be a vector', diag.type) 1802 1803 return gof.Apply(self, [diag], 1804 [SparseType(dtype=diag.dtype, format='csc')()]) 1805 1806 def perform(self, node, inputs, outputs): 1807 (z,) = outputs 1808 diag = inputs[0] 1809 1810 N = len(diag) 1811 data = diag[:N] 1812 indices = list(range(N)) 1813 indptr = list(range(N + 1)) 1814 tup = (data, indices, indptr) 1815 1816 z[0] = scipy.sparse.csc_matrix(tup, copy=True) 1817 1818 def grad(self, inputs, gout): 1819 (gz,) = gout 1820 return [diag(gz)] 1821 1822 def infer_shape(self, nodes, shapes): 1823 return [(shapes[0][0], shapes[0][0])] 1824 1825square_diagonal = SquareDiagonal() 1826""" 1827Return a square sparse (csc) matrix whose diagonal is given by the dense vector 1828argument. 1829 1830Parameters 1831---------- 1832x 1833 Dense vector for the diagonal. 1834 1835Returns 1836------- 1837sparse matrix 1838 A sparse matrix having `x` as diagonal. 1839 1840Notes 1841----- 1842The grad implemented is regular, i.e. not structured. 1843 1844""" 1845 1846 1847class EnsureSortedIndices(gof.op.Op): 1848 # See doc in instance of this Op or function after this class definition. 1849 __props__ = ("inplace",) 1850 1851 def __init__(self, inplace): 1852 self.inplace = inplace 1853 if self.inplace: 1854 self.view_map = {0: [0]} 1855 1856 def make_node(self, x): 1857 x = as_sparse_variable(x) 1858 assert x.format in ["csr", "csc"] 1859 return gof.Apply(self, [x], [x.type()]) 1860 1861 def perform(self, node, inputs, outputs): 1862 (x,) = inputs 1863 (z,) = outputs 1864 if self.inplace: 1865 z[0] = x.sort_indices() 1866 else: 1867 z[0] = x.sorted_indices() 1868 1869 def grad(self, inputs, output_grad): 1870 return [output_grad[0]] 1871 1872 def infer_shape(self, node, i0_shapes): 1873 return i0_shapes 1874 1875 def __str__(self): 1876 if self.inplace: 1877 return self.__class__.__name__ + "{inplace}" 1878 else: 1879 return self.__class__.__name__ + "{no_inplace}" 1880ensure_sorted_indices = EnsureSortedIndices(inplace=False) 1881""" 1882Re-sort indices of a sparse matrix. 1883 1884CSR column indices are not necessarily sorted. Likewise 1885for CSC row indices. Use `ensure_sorted_indices` when sorted 1886indices are required (e.g. when passing data to other 1887libraries). 1888 1889Parameters 1890---------- 1891x 1892 A sparse matrix. 1893 1894Returns 1895------- 1896sparse matrix 1897 The same as `x` with indices sorted. 1898 1899Notes 1900----- 1901The grad implemented is regular, i.e. not structured. 1902 1903""" 1904 1905 1906def clean(x): 1907 """ 1908 Remove explicit zeros from a sparse matrix, and re-sort indices. 1909 1910 CSR column indices are not necessarily sorted. Likewise 1911 for CSC row indices. Use `clean` when sorted 1912 indices are required (e.g. when passing data to other 1913 libraries) and to ensure there are no zeros in the data. 1914 1915 Parameters 1916 ---------- 1917 x 1918 A sparse matrix. 1919 1920 Returns 1921 ------- 1922 A sparse matrix 1923 The same as `x` with indices sorted and zeros 1924 removed. 1925 1926 Notes 1927 ----- 1928 The grad implemented is regular, i.e. not structured. 1929 1930 """ 1931 return ensure_sorted_indices(remove0(x)) 1932 1933 1934class AddSS(gof.op.Op): 1935 # add(sparse, sparse). 1936 # see the doc of add() for more detail. 1937 __props__ = () 1938 1939 def make_node(self, x, y): 1940 x, y = map(as_sparse_variable, [x, y]) 1941 assert x.format in ["csr", "csc"] 1942 assert y.format in ["csr", "csc"] 1943 out_dtype = scalar.upcast(x.type.dtype, y.type.dtype) 1944 return gof.Apply(self, 1945 [x, y], 1946 [SparseType(dtype=out_dtype, 1947 format=x.type.format)()]) 1948 1949 def perform(self, node, inputs, outputs): 1950 (x, y) = inputs 1951 (out,) = outputs 1952 assert _is_sparse(x) and _is_sparse(y) 1953 assert x.shape == y.shape 1954 out[0] = x + y 1955 1956 def grad(self, inputs, gout): 1957 (x, y) = inputs 1958 (gz,) = gout 1959 assert _is_sparse_variable(x) and _is_sparse_variable(y) 1960 assert _is_sparse_variable(gz) 1961 return gz, gz 1962 1963 def infer_shape(self, node, shapes): 1964 return [shapes[0]] 1965 1966add_s_s = AddSS() 1967 1968 1969class AddSSData(gof.op.Op): 1970 # See doc in instance of this Op or function after this class definition. 1971 __props__ = () 1972 1973 def make_node(self, x, y): 1974 x, y = map(as_sparse_variable, [x, y]) 1975 assert x.format in ["csr", "csc"] 1976 assert y.format in ["csr", "csc"] 1977 if x.type.dtype != y.type.dtype: 1978 raise NotImplementedError() 1979 if x.type.format != y.type.format: 1980 raise NotImplementedError() 1981 return gof.Apply(self, 1982 [x, y], 1983 [SparseType(dtype=x.type.dtype, 1984 format=x.type.format)()]) 1985 1986 def perform(self, node, inputs, outputs): 1987 (x, y) = inputs 1988 (out,) = outputs 1989 assert _is_sparse(x) and _is_sparse(y) 1990 assert x.shape == y.shape 1991 assert x.data.shape == y.data.shape 1992 out[0] = x.copy() 1993 out[0].data += y.data 1994 1995 def grad(self, inputs, gout): 1996 (gz,) = gout 1997 is_continuous = [(i.dtype in continuous_dtypes) 1998 for i in inputs] 1999 derivative = {True: gz, False: None} 2000 return [derivative[b] for b in is_continuous] 2001 2002 def infer_shape(self, node, ins_shapes): 2003 return [ins_shapes[0]] 2004 2005add_s_s_data = AddSSData() 2006""" 2007Add two sparse matrices assuming they have the same sparsity pattern. 2008 2009Parameters 2010---------- 2011x 2012 Sparse matrix. 2013y 2014 Sparse matrix. 2015 2016Returns 2017------- 2018A sparse matrix 2019 The sum of the two sparse matrices element wise. 2020 2021Notes 2022----- 2023`x` and `y` are assumed to have the same sparsity pattern. 2024 2025The grad implemented is structured. 2026 2027""" 2028 2029 2030class AddSD(gof.op.Op): 2031 # add(sparse, sparse). 2032 # see the doc of add() for more detail. 2033 __props__ = () 2034 2035 def make_node(self, x, y): 2036 x, y = as_sparse_variable(x), tensor.as_tensor_variable(y) 2037 assert x.format in ["csr", "csc"] 2038 out_dtype = scalar.upcast(x.type.dtype, y.type.dtype) 2039 2040 # The magic number two here arises because L{scipy.sparse} 2041 # objects must be matrices (have dimension 2) 2042 assert y.type.ndim == 2 2043 return gof.Apply(self, 2044 [x, y], 2045 [tensor.TensorType(dtype=out_dtype, 2046 broadcastable=y.type.broadcastable 2047 )()]) 2048 2049 def perform(self, node, inputs, outputs): 2050 (x, y) = inputs 2051 (out,) = outputs 2052 assert _is_dense(y) 2053 2054 # The asarray is needed as in some case, this return a 2055 # numpy.matrixlib.defmatrix.matrix object and not an ndarray. 2056 out[0] = theano._asarray(x + y, dtype=node.outputs[0].type.dtype) 2057 2058 def grad(self, inputs, gout): 2059 (x, y) = inputs 2060 (gz,) = gout 2061 assert _is_sparse_variable(x) and _is_dense_variable(y) 2062 assert _is_dense_variable(gz) 2063 return sp_ones_like(x) * gz, gz 2064 2065 def infer_shape(self, node, shapes): 2066 return [shapes[1]] 2067 2068add_s_d = AddSD() 2069 2070 2071class StructuredAddSV(gof.op.Op): 2072 2073 __props__ = () 2074 2075 def make_node(self, x, y): 2076 x = as_sparse_variable(x) 2077 assert x.format in ["csr", "csc"] 2078 y = tensor.as_tensor_variable(y) 2079 2080 assert y.type.ndim == 1 2081 2082 if x.type.dtype != y.type.dtype: 2083 raise NotImplementedError() 2084 return gof.Apply(self, 2085 [x, y], 2086 [SparseType(dtype=x.type.dtype, 2087 format=x.type.format)()]) 2088 2089 def perform(self, node, inputs, outputs): 2090 (x, y) = inputs 2091 (out,) = outputs 2092 assert _is_sparse(x) and not _is_sparse(y) 2093 assert x.shape[1] == y.shape[0] 2094 out[0] = x.__class__(x + (x.toarray() != 0) * y) 2095 2096 def grad(self, inputs, gout): 2097 (x, y) = inputs 2098 (gz,) = gout 2099 assert _is_sparse_variable(x) and not _is_sparse_variable(y) 2100 assert _is_sparse_variable(gz) 2101 return gz, sp_sum(gz, axis=0, sparse_grad=True) 2102 2103 def infer_shape(self, node, ins_shapes): 2104 return [ins_shapes[0]] 2105 2106structured_add_s_v = StructuredAddSV() 2107""" 2108Structured addition of a sparse matrix and a dense vector. 2109The elements of the vector are only added to the corresponding 2110non-zero elements of the sparse matrix. Therefore, this operation 2111outputs another sparse matrix. 2112 2113Parameters 2114---------- 2115x 2116 Sparse matrix. 2117y 2118 Tensor type vector. 2119 2120Returns 2121------- 2122A sparse matrix 2123 A sparse matrix containing the addition of the vector to 2124 the data of the sparse matrix. 2125 2126Notes 2127----- 2128The grad implemented is structured since the op is structured. 2129 2130""" 2131 2132 2133def add(x, y): 2134 """ 2135 Add two matrices, at least one of which is sparse. 2136 2137 This method will provide the right op according 2138 to the inputs. 2139 2140 Parameters 2141 ---------- 2142 x 2143 A matrix variable. 2144 y 2145 A matrix variable. 2146 2147 Returns 2148 ------- 2149 A sparse matrix 2150 `x` + `y` 2151 2152 Notes 2153 ----- 2154 At least one of `x` and `y` must be a sparse matrix. 2155 2156 The grad will be structured only when one of the variable will be a dense 2157 matrix. 2158 2159 """ 2160 2161 if hasattr(x, 'getnnz'): 2162 x = as_sparse_variable(x) 2163 if hasattr(y, 'getnnz'): 2164 y = as_sparse_variable(y) 2165 if not isinstance(x, theano.Variable): 2166 x = theano.tensor.as_tensor_variable(x) 2167 if not isinstance(y, theano.Variable): 2168 y = theano.tensor.as_tensor_variable(y) 2169 2170 x_is_sparse_variable = _is_sparse_variable(x) 2171 y_is_sparse_variable = _is_sparse_variable(y) 2172 2173 assert x_is_sparse_variable or y_is_sparse_variable 2174 if x_is_sparse_variable and y_is_sparse_variable: 2175 return add_s_s(x, y) 2176 elif x_is_sparse_variable and not y_is_sparse_variable: 2177 return add_s_d(x, y) 2178 elif y_is_sparse_variable and not x_is_sparse_variable: 2179 return add_s_d(y, x) 2180 else: 2181 raise NotImplementedError() 2182 2183 2184def sub(x, y): 2185 """ 2186 Subtract two matrices, at least one of which is sparse. 2187 2188 This method will provide the right op according 2189 to the inputs. 2190 2191 Parameters 2192 ---------- 2193 x 2194 A matrix variable. 2195 y 2196 A matrix variable. 2197 2198 Returns 2199 ------- 2200 A sparse matrix 2201 `x` - `y` 2202 2203 Notes 2204 ----- 2205 At least one of `x` and `y` must be a sparse matrix. 2206 2207 The grad will be structured only when one of the variable will be a dense 2208 matrix. 2209 2210 """ 2211 return x + (-y) 2212 2213 2214class MulSS(gof.op.Op): 2215 # mul(sparse, sparse) 2216 # See the doc of mul() for more detail 2217 __props__ = () 2218 2219 def make_node(self, x, y): 2220 x, y = as_sparse_variable(x), as_sparse_variable(y) 2221 assert x.format in ["csr", "csc"] 2222 assert y.format in ["csr", "csc"] 2223 out_dtype = scalar.upcast(x.type.dtype, y.type.dtype) 2224 return gof.Apply(self, 2225 [x, y], 2226 [SparseType(dtype=out_dtype, 2227 format=x.type.format)()]) 2228 2229 def perform(self, node, inputs, outputs): 2230 (x, y) = inputs 2231 (out,) = outputs 2232 assert _is_sparse(x) and _is_sparse(y) 2233 assert len(x.shape) == 2 2234 assert y.shape == x.shape 2235 # This calls the element-wise multiple 2236 # x * y calls dot... 2237 out[0] = x.multiply(y) 2238 2239 def grad(self, inputs, gout): 2240 (x, y) = inputs 2241 (gz,) = gout 2242 return y * gz, x * gz 2243 2244 def infer_shape(self, node, shapes): 2245 return [shapes[0]] 2246 2247mul_s_s = MulSS() 2248 2249 2250class MulSD(gof.op.Op): 2251 # mul(sparse, dense) 2252 # See the doc of mul() for more detail 2253 __props__ = () 2254 2255 def make_node(self, x, y): 2256 x, y = as_sparse_variable(x), tensor.as_tensor_variable(y) 2257 2258 assert x.format in ["csr", "csc"] 2259 2260 # upcast the tensor. Is the cast of sparse done implemented? 2261 dtype = scalar.upcast(x.type.dtype, y.type.dtype) 2262 2263 # The magic number two here arises because L{scipy.sparse} 2264 # objects must be matrices (have dimension 2) 2265 # Broadcasting of the sparse matrix is not supported. 2266 # We support nd == 0 used by grad of SpSum() 2267 assert y.type.ndim in [0, 2] 2268 out = SparseType(dtype=dtype, 2269 format=x.type.format)() 2270 return gof.Apply(self, [x, y], [out]) 2271 2272 def perform(self, node, inputs, outputs): 2273 (x, y) = inputs 2274 (out,) = outputs 2275 assert _is_sparse(x) and _is_dense(y) 2276 if len(y.shape) == 0: 2277 out_dtype = node.outputs[0].dtype 2278 if x.dtype == out_dtype: 2279 z = x.copy() 2280 else: 2281 z = x.astype(out_dtype) 2282 out[0] = z 2283 out[0].data *= y 2284 elif len(y.shape) == 1: 2285 raise NotImplementedError() # RowScale / ColScale 2286 elif len(y.shape) == 2: 2287 # if we have enough memory to fit y, maybe we can fit x.asarray() 2288 # too? 2289 # TODO: change runtime from O(M*N) to O(nonzeros) 2290 M, N = x.shape 2291 assert x.shape == y.shape 2292 out_dtype = node.outputs[0].dtype 2293 2294 if x.format == 'csc': 2295 indices = x.indices 2296 indptr = x.indptr 2297 if x.dtype == out_dtype: 2298 z = x.copy() 2299 else: 2300 z = x.astype(out_dtype) 2301 z_data = z.data 2302 2303 for j in xrange(0, N): 2304 for i_idx in xrange(indptr[j], indptr[j + 1]): 2305 i = indices[i_idx] 2306 z_data[i_idx] *= y[i, j] 2307 out[0] = z 2308 elif x.format == 'csr': 2309 indices = x.indices 2310 indptr = x.indptr 2311 if x.dtype == out_dtype: 2312 z = x.copy() 2313 else: 2314 z = x.astype(out_dtype) 2315 z_data = z.data 2316 2317 for i in xrange(0, M): 2318 for j_idx in xrange(indptr[i], indptr[i + 1]): 2319 j = indices[j_idx] 2320 z_data[j_idx] *= y[i, j] 2321 out[0] = z 2322 else: 2323 print(( 2324 "WARNING: crappy implementation of MulSD" 2325 ), x.format, file=sys.stderr) 2326 out[0] = type(x)(x.toarray() * y) 2327 2328 def grad(self, inputs, gout): 2329 (x, y) = inputs 2330 (gz,) = gout 2331 assert _is_sparse_variable(x) and _is_dense_variable(y) 2332 assert _is_sparse_variable(gz) 2333 return y * gz, dense_from_sparse(x * gz) 2334 2335 def infer_shape(self, node, shapes): 2336 return [shapes[0]] 2337mul_s_d = MulSD() 2338 2339 2340class MulSV(gof.op.Op): 2341 2342 __props__ = () 2343 2344 def make_node(self, x, y): 2345 x = as_sparse_variable(x) 2346 assert x.format in ["csr", "csc"] 2347 y = tensor.as_tensor_variable(y) 2348 2349 assert y.type.ndim == 1 2350 2351 if x.type.dtype != y.type.dtype: 2352 raise NotImplementedError( 2353 "MulSV not implemented for differing dtypes." 2354 "Got %s and %s." % (str(x.type.dtype), str(y.type.dtype))) 2355 return gof.Apply(self, 2356 [x, y], 2357 [SparseType(dtype=x.type.dtype, 2358 format=x.type.format)()]) 2359 2360 def perform(self, node, inputs, outputs): 2361 (x, y) = inputs 2362 (out,) = outputs 2363 assert _is_sparse(x) and not _is_sparse(y) 2364 assert x.shape[1] == y.shape[0] 2365 out[0] = x.__class__(x.toarray() * y) 2366 2367 def grad(self, inputs, gout): 2368 (x, y) = inputs 2369 (gz,) = gout 2370 assert _is_sparse_variable(x) and _is_dense_variable(y) 2371 assert _is_sparse_variable(gz) 2372 2373 # mul_s_v is not implemented if the types vary 2374 2375 if gz.dtype == 'float64' and y.dtype == 'float32': 2376 y = y.astype('float64') 2377 2378 if gz.dtype == 'float32' and y.dtype == 'float64': 2379 gz = gz.astype('float64') 2380 2381 return mul_s_v(gz, y), sp_sum(x * gz, axis=0, sparse_grad=True) 2382 2383 def infer_shape(self, node, ins_shapes): 2384 return [ins_shapes[0]] 2385 2386mul_s_v = MulSV() 2387""" 2388Multiplication of sparse matrix by a broadcasted dense vector element wise. 2389 2390Parameters 2391---------- 2392x 2393 Sparse matrix to multiply. 2394y 2395 Tensor broadcastable vector. 2396 2397Returns 2398------- 2399A sparse matrix 2400 The product x * y element wise. 2401 2402Notes 2403----- 2404The grad implemented is regular, i.e. not structured. 2405 2406""" 2407 2408 2409def mul(x, y): 2410 """ 2411 Multiply elementwise two matrices, at least one of which is sparse. 2412 2413 This method will provide the right op according to the inputs. 2414 2415 Parameters 2416 ---------- 2417 x 2418 A matrix variable. 2419 y 2420 A matrix variable. 2421 2422 Returns 2423 ------- 2424 A sparse matrix 2425 `x` * `y` 2426 2427 Notes 2428 ----- 2429 At least one of `x` and `y` must be a sparse matrix. 2430 The grad is regular, i.e. not structured. 2431 2432 """ 2433 2434 x = as_sparse_or_tensor_variable(x) 2435 y = as_sparse_or_tensor_variable(y) 2436 2437 x_is_sparse_variable = _is_sparse_variable(x) 2438 y_is_sparse_variable = _is_sparse_variable(y) 2439 2440 assert x_is_sparse_variable or y_is_sparse_variable 2441 if x_is_sparse_variable and y_is_sparse_variable: 2442 2443 # mul_s_s is not implemented if the types differ 2444 if y.dtype == 'float64' and x.dtype == 'float32': 2445 x = x.astype('float64') 2446 2447 return mul_s_s(x, y) 2448 elif x_is_sparse_variable and not y_is_sparse_variable: 2449 2450 # mul is unimplemented if the dtypes differ 2451 if y.dtype == 'float64' and x.dtype == 'float32': 2452 x = x.astype('float64') 2453 2454 return mul_s_d(x, y) 2455 elif y_is_sparse_variable and not x_is_sparse_variable: 2456 return mul_s_d(y, x) 2457 else: 2458 raise NotImplementedError() 2459 2460 2461class __ComparisonOpSS(gof.op.Op): 2462 """ 2463 Used as a superclass for all comparisons between two sparses matrices. 2464 2465 Parameters 2466 ---------- 2467 x 2468 First compared sparse matrix. 2469 y 2470 Second compared sparse matrix 2471 2472 Returns 2473 ------- 2474 object 2475 Comparison(x,y) 2476 2477 """ 2478 2479 __props__ = () 2480 2481 # Function to override 2482 def comparison(self, x, y): 2483 raise NotImplementedError() 2484 2485 def make_node(self, x, y): 2486 x = as_sparse_variable(x) 2487 y = as_sparse_variable(y) 2488 2489 if x.type.format != y.type.format: 2490 raise NotImplementedError() 2491 return gof.Apply(self, 2492 [x, y], 2493 [SparseType(dtype='uint8', 2494 format=x.type.format)()]) 2495 2496 def perform(self, node, inputs, outputs): 2497 (x, y) = inputs 2498 (out,) = outputs 2499 assert _is_sparse(x) and _is_sparse(y) 2500 assert x.shape == y.shape 2501 out[0] = self.comparison(x, y).astype('uint8') 2502 2503 def infer_shape(self, node, ins_shapes): 2504 return [ins_shapes[0]] 2505 2506 2507class __ComparisonOpSD(gof.op.Op): 2508 """ 2509 Used as a superclass for all comparisons between sparse and dense matrix. 2510 2511 Parameters 2512 ---------- 2513 x 2514 Sparse matrix. 2515 y 2516 Dense matrix. 2517 2518 Returns 2519 ------- 2520 object 2521 Comparison(x,y) 2522 2523 """ 2524 2525 __props__ = () 2526 2527 # Function to override 2528 def comparison(self, x, y): 2529 raise NotImplementedError() 2530 2531 def make_node(self, x, y): 2532 x, y = as_sparse_variable(x), tensor.as_tensor_variable(y) 2533 2534 assert y.type.ndim == 2 2535 out = tensor.TensorType(dtype='uint8', broadcastable=(False, False))() 2536 return gof.Apply(self, 2537 [x, y], 2538 [out]) 2539 2540 def perform(self, node, inputs, outputs): 2541 (x, y) = inputs 2542 (out,) = outputs 2543 assert _is_sparse(x) 2544 assert x.shape == y.shape 2545 assert _is_dense(y) 2546 o = self.comparison(x, y).astype('uint8') 2547 o = np.asarray(o) 2548 out[0] = o 2549 2550 def infer_shape(self, node, ins_shapes): 2551 return [ins_shapes[0]] 2552 2553 2554def __ComparisonSwitch(SS, SD, DS): 2555 """ 2556 2557 Parameters 2558 ---------- 2559 SS 2560 Function to apply between two sparses matrices. 2561 SD 2562 Function to apply between a sparse and a dense matrix. 2563 DS 2564 Function to apply between a dense and a sparse matrix. 2565 2566 Returns 2567 ------- 2568 function 2569 Switch function taking two matrices as input. 2570 2571 Notes 2572 ----- 2573 At least one of `x` and `y` must be a sparse matrix. 2574 2575 DS swap input as a dense matrix cannot be a left operand. 2576 2577 """ 2578 2579 def helper(x, y): 2580 2581 scipy_ver = [int(n) for n in scipy.__version__.split('.')[:2]] 2582 2583 assert scipy_ver >= [0, 13] 2584 2585 if hasattr(x, 'getnnz'): 2586 x = as_sparse_variable(x) 2587 if hasattr(y, 'getnnz'): 2588 y = as_sparse_variable(y) 2589 if not isinstance(x, theano.Variable): 2590 x = theano.tensor.as_tensor_variable(x) 2591 if not isinstance(y, theano.Variable): 2592 y = theano.tensor.as_tensor_variable(y) 2593 2594 x_is_sparse_variable = _is_sparse_variable(x) 2595 y_is_sparse_variable = _is_sparse_variable(y) 2596 2597 assert x_is_sparse_variable or y_is_sparse_variable 2598 if x_is_sparse_variable and y_is_sparse_variable: 2599 return SS(x, y) 2600 elif x_is_sparse_variable and not y_is_sparse_variable: 2601 return SD(x, y) 2602 elif y_is_sparse_variable and not x_is_sparse_variable: 2603 return DS(y, x) 2604 else: 2605 raise NotImplementedError() 2606 2607 return helper 2608 2609 2610class EqualSS(__ComparisonOpSS): 2611 def comparison(self, x, y): 2612 return x == y 2613 2614 2615equal_s_s = EqualSS() 2616 2617 2618class EqualSD(__ComparisonOpSD): 2619 def comparison(self, x, y): 2620 return x == y 2621 2622equal_s_d = EqualSD() 2623 2624 2625class NotEqualSS(__ComparisonOpSS): 2626 def comparison(self, x, y): 2627 return x != y 2628 2629not_equal_s_s = NotEqualSS() 2630 2631 2632class NotEqualSD(__ComparisonOpSD): 2633 def comparison(self, x, y): 2634 return x != y 2635 2636not_equal_s_d = NotEqualSD() 2637 2638 2639class LessThanSS(__ComparisonOpSS): 2640 def comparison(self, x, y): 2641 return x < y 2642 2643less_than_s_s = LessThanSS() 2644 2645 2646class LessThanSD(__ComparisonOpSD): 2647 def comparison(self, x, y): 2648 return x < y 2649 2650less_than_s_d = LessThanSD() 2651 2652 2653class GreaterThanSS(__ComparisonOpSS): 2654 def comparison(self, x, y): 2655 return x > y 2656 2657greater_than_s_s = GreaterThanSS() 2658 2659 2660class GreaterThanSD(__ComparisonOpSD): 2661 def comparison(self, x, y): 2662 return x > y 2663 2664greater_than_s_d = GreaterThanSD() 2665 2666 2667class LessEqualSS(__ComparisonOpSS): 2668 def comparison(self, x, y): 2669 return x <= y 2670 2671less_equal_s_s = LessEqualSS() 2672 2673 2674class LessEqualSD(__ComparisonOpSD): 2675 def comparison(self, x, y): 2676 return x <= y 2677 2678less_equal_s_d = LessEqualSD() 2679 2680 2681class GreaterEqualSS(__ComparisonOpSS): 2682 def comparison(self, x, y): 2683 return x >= y 2684 2685greater_equal_s_s = GreaterEqualSS() 2686 2687 2688class GreaterEqualSD(__ComparisonOpSD): 2689 def comparison(self, x, y): 2690 return x >= y 2691 2692greater_equal_s_d = GreaterEqualSD() 2693 2694 2695eq = __ComparisonSwitch(equal_s_s, equal_s_d, equal_s_d) 2696""" 2697Parameters 2698---------- 2699x 2700 A matrix variable. 2701y 2702 A matrix variable. 2703 2704Returns 2705------- 2706matrix variable 2707 `x` == `y` 2708 2709Notes 2710----- 2711At least one of `x` and `y` must be a sparse matrix. 2712 2713""" 2714 2715 2716neq = __ComparisonSwitch(not_equal_s_s, not_equal_s_d, not_equal_s_d) 2717""" 2718Parameters 2719---------- 2720x 2721 A matrix variable. 2722y 2723 A matrix variable. 2724 2725Returns 2726------- 2727matrix variable 2728 `x` != `y` 2729 2730Notes 2731----- 2732At least one of `x` and `y` must be a sparse matrix. 2733 2734""" 2735 2736 2737lt = __ComparisonSwitch(less_than_s_s, less_than_s_d, greater_than_s_d) 2738""" 2739Parameters 2740---------- 2741x 2742 A matrix variable. 2743y 2744 A matrix variable. 2745 2746Returns 2747------- 2748matrix variable 2749 `x` < `y` 2750 2751Notes 2752----- 2753At least one of `x` and `y` must be a sparse matrix. 2754 2755""" 2756 2757 2758gt = __ComparisonSwitch(greater_than_s_s, greater_than_s_d, less_than_s_d) 2759""" 2760Parameters 2761---------- 2762x 2763 A matrix variable. 2764y 2765 A matrix variable. 2766 2767Returns 2768------- 2769matrix variable 2770 `x` > `y` 2771 2772Notes 2773----- 2774At least one of `x` and `y` must be a sparse matrix. 2775 2776""" 2777 2778le = __ComparisonSwitch(less_equal_s_s, less_equal_s_d, greater_equal_s_d) 2779""" 2780Parameters 2781---------- 2782x 2783 A matrix variable. 2784y 2785 A matrix variable. 2786 2787Returns 2788------- 2789matrix variable 2790 `x` <= `y` 2791 2792Notes 2793----- 2794At least one of `x` and `y` must be a sparse matrix. 2795 2796""" 2797 2798ge = __ComparisonSwitch(greater_equal_s_s, greater_equal_s_d, 2799 less_equal_s_d) 2800""" 2801Parameters 2802---------- 2803x 2804 A matrix variable. 2805y 2806 A matrix variable. 2807 2808Returns 2809------- 2810matrix variable 2811 `x` >= `y` 2812 2813Notes 2814----- 2815At least one of `x` and `y` must be a sparse matrix. 2816 2817""" 2818 2819 2820class HStack(gof.op.Op): 2821 # See doc in instance of this Op or function after this class definition. 2822 __props__ = ("format", "dtype") 2823 2824 def __init__(self, format=None, dtype=None): 2825 if format is None: 2826 self.format = 'csc' 2827 else: 2828 self.format = format 2829 2830 if dtype is None: 2831 raise ValueError('The output dtype must be specified.') 2832 self.dtype = dtype 2833 2834 def make_node(self, *mat): 2835 if not mat: 2836 raise ValueError('Cannot join an empty list of sparses.') 2837 var = [as_sparse_variable(x) for x in mat] 2838 2839 for x in var: 2840 assert x.format in ["csr", "csc"] 2841 2842 return gof.Apply(self, 2843 var, 2844 [SparseType(dtype=self.dtype, 2845 format=self.format)()]) 2846 2847 def perform(self, node, block, outputs): 2848 (out,) = outputs 2849 for b in block: 2850 assert _is_sparse(b) 2851 out[0] = scipy.sparse.hstack(block, format=self.format, 2852 dtype=self.dtype) 2853 # Some version of scipy (at least 0.14.0.dev-c4314b0) 2854 # Do not cast to the wanted dtype. 2855 if out[0].dtype != self.dtype: 2856 out[0] = out[0].astype(self.dtype) 2857 2858 def grad(self, inputs, gout): 2859 (gz,) = gout 2860 is_continuous = [(inputs[i].dtype in tensor.continuous_dtypes) 2861 for i in range(len(inputs))] 2862 2863 if _is_sparse_variable(gz): 2864 gz = dense_from_sparse(gz) 2865 2866 split = tensor.Split(len(inputs))(gz, 1, 2867 tensor.stack( 2868 [x.shape[1] 2869 for x in inputs])) 2870 if not isinstance(split, list): 2871 split = [split] 2872 2873 derivative = [SparseFromDense(self.format)(s) for s in split] 2874 2875 def choose(continuous, derivative): 2876 if continuous: 2877 return derivative 2878 else: 2879 return None 2880 return [choose(c, d) for c, d in zip(is_continuous, derivative)] 2881 2882 def infer_shape(self, node, ins_shapes): 2883 def _get(l): 2884 return l[1] 2885 d = sum(map(_get, ins_shapes)) 2886 return [(ins_shapes[0][0], d)] 2887 2888 def __str__(self): 2889 return "%s(%s,%s)" % (self.__class__.__name__, self.format, self.dtype) 2890 2891 2892def hstack(blocks, format=None, dtype=None): 2893 """ 2894 Stack sparse matrices horizontally (column wise). 2895 2896 This wrap the method hstack from scipy. 2897 2898 Parameters 2899 ---------- 2900 blocks 2901 List of sparse array of compatible shape. 2902 format 2903 String representing the output format. Default is csc. 2904 dtype 2905 Output dtype. 2906 2907 Returns 2908 ------- 2909 array 2910 The concatenation of the sparse array column wise. 2911 2912 Notes 2913 ----- 2914 The number of line of the sparse matrix must agree. 2915 2916 The grad implemented is regular, i.e. not structured. 2917 2918 """ 2919 2920 blocks = [as_sparse_variable(i) for i in blocks] 2921 if dtype is None: 2922 dtype = theano.scalar.upcast(*[i.dtype for i in blocks]) 2923 return HStack(format=format, dtype=dtype)(*blocks) 2924 2925 2926class VStack(HStack): 2927 # See doc in instance of this Op or function after this class definition. 2928 def perform(self, node, block, outputs): 2929 (out,) = outputs 2930 for b in block: 2931 assert _is_sparse(b) 2932 out[0] = scipy.sparse.vstack(block, format=self.format, 2933 dtype=self.dtype) 2934 # Some version of scipy (at least 0.14.0.dev-c4314b0) 2935 # Do not cast to the wanted dtype. 2936 if out[0].dtype != self.dtype: 2937 out[0] = out[0].astype(self.dtype) 2938 2939 def grad(self, inputs, gout): 2940 (gz,) = gout 2941 is_continuous = [(inputs[i].dtype in tensor.continuous_dtypes) 2942 for i in range(len(inputs))] 2943 2944 if _is_sparse_variable(gz): 2945 gz = dense_from_sparse(gz) 2946 2947 split = tensor.Split(len(inputs))(gz, 0, 2948 tensor.stack( 2949 [x.shape[0] 2950 for x in inputs])) 2951 if not isinstance(split, list): 2952 split = [split] 2953 2954 derivative = [SparseFromDense(self.format)(s) for s in split] 2955 2956 def choose(continuous, derivative): 2957 if continuous: 2958 return derivative 2959 else: 2960 return None 2961 return [choose(c, d) for c, d in zip(is_continuous, derivative)] 2962 2963 def infer_shape(self, node, ins_shapes): 2964 def _get(l): 2965 return l[0] 2966 d = sum(map(_get, ins_shapes)) 2967 return [(d, ins_shapes[0][1])] 2968 2969 2970def vstack(blocks, format=None, dtype=None): 2971 """ 2972 Stack sparse matrices vertically (row wise). 2973 2974 This wrap the method vstack from scipy. 2975 2976 Parameters 2977 ---------- 2978 blocks 2979 List of sparse array of compatible shape. 2980 format 2981 String representing the output format. Default is csc. 2982 dtype 2983 Output dtype. 2984 2985 Returns 2986 ------- 2987 array 2988 The concatenation of the sparse array row wise. 2989 2990 Notes 2991 ----- 2992 The number of column of the sparse matrix must agree. 2993 2994 The grad implemented is regular, i.e. not structured. 2995 2996 """ 2997 2998 blocks = [as_sparse_variable(i) for i in blocks] 2999 if dtype is None: 3000 dtype = theano.scalar.upcast(*[i.dtype for i in blocks]) 3001 return VStack(format=format, dtype=dtype)(*blocks) 3002 3003 3004class Remove0(gof.Op): 3005 # See doc in instance of this Op or a function after the class definition. 3006 __props__ = ("inplace",) 3007 3008 def __init__(self, inplace=False): 3009 self.inplace = inplace 3010 if self.inplace: 3011 self.destroy_map = {0: [0]} 3012 3013 def __str__(self): 3014 l = [] 3015 if self.inplace: 3016 l.append('inplace') 3017 return self.__class__.__name__ + '{%s}' % ', '.join(l) 3018 3019 def make_node(self, x): 3020 x = as_sparse_variable(x) 3021 assert x.format in ["csr", "csc"] 3022 return gof.Apply(self, [x], [x.type()]) 3023 3024 def perform(self, node, inputs, outputs): 3025 (x,) = inputs 3026 (z,) = outputs 3027 if self.inplace: 3028 c = x 3029 else: 3030 c = x.copy() 3031 c.eliminate_zeros() 3032 z[0] = c 3033 3034 def grad(self, inputs, gout): 3035 (x,) = inputs 3036 (gz,) = gout 3037 return [gz] 3038 3039 def infer_shape(self, node, i0_shapes): 3040 return i0_shapes 3041remove0 = Remove0() 3042""" 3043Remove explicit zeros from a sparse matrix. 3044 3045Parameters 3046---------- 3047x 3048 Sparse matrix. 3049 3050Returns 3051------- 3052sparse matrix 3053 Exactly `x` but with a data attribute exempt of zeros. 3054 3055Notes 3056----- 3057The grad implemented is regular, i.e. not structured. 3058 3059""" 3060 3061 3062# Structured monoid 3063def structured_monoid(tensor_op): 3064 # Generic operation to perform many kinds of monoid element-wise 3065 # operations on the non-zeros of a sparse matrix. 3066 3067 # The first parameter must always be a sparse matrix. The other parameters 3068 # must be scalars which will be passed as argument to the tensor_op. 3069 3070 def decorator(f): 3071 def wrapper(*args): 3072 x = as_sparse_variable(args[0]) 3073 assert x.format in ["csr", "csc"] 3074 3075 xs = [scalar.as_scalar(arg) for arg in args[1:]] 3076 3077 data, ind, ptr, shape = csm_properties(x) 3078 3079 data = tensor_op(data, *xs) 3080 3081 return CSM(x.format)(data, ind, ptr, shape) 3082 wrapper.__name__ = str(tensor_op.scalar_op) 3083 return wrapper 3084 return decorator 3085 3086 3087@structured_monoid(tensor.nnet.sigmoid) 3088def structured_sigmoid(x): 3089 """ 3090 Structured elemwise sigmoid. 3091 3092 """ 3093 # see decorator for function body 3094 3095 3096@structured_monoid(tensor.exp) 3097def structured_exp(x): 3098 """ 3099 Structured elemwise exponential. 3100 3101 """ 3102 # see decorator for function body 3103 3104 3105@structured_monoid(tensor.log) 3106def structured_log(x): 3107 """ 3108 Structured elemwise logarithm. 3109 3110 """ 3111 # see decorator for function body 3112 3113 3114@structured_monoid(tensor.pow) 3115def structured_pow(x, y): 3116 """ 3117 Structured elemwise power of sparse matrix x by scalar y. 3118 3119 """ 3120 # see decorator for function body 3121 3122 3123@structured_monoid(tensor.minimum) 3124def structured_minimum(x, y): 3125 """ 3126 Structured elemwise minimum of sparse matrix x by scalar y. 3127 3128 """ 3129 # see decorator for function body 3130 3131 3132@structured_monoid(tensor.maximum) 3133def structured_maximum(x, y): 3134 """ 3135 Structured elemwise maximum of sparse matrix x by scalar y. 3136 3137 """ 3138 # see decorator for function body 3139 3140 3141@structured_monoid(tensor.add) 3142def structured_add(x): 3143 """ 3144 Structured addition of sparse matrix x and scalar y. 3145 3146 """ 3147 # see decorator for function body 3148 3149 3150# Sparse operation (map 0 to 0) 3151@structured_monoid(tensor.sin) 3152def sin(x): 3153 """ 3154 Elemwise sinus of `x`. 3155 3156 """ 3157 # see decorator for function body 3158 3159 3160@structured_monoid(tensor.tan) 3161def tan(x): 3162 """ 3163 Elemwise tan of `x`. 3164 3165 """ 3166 # see decorator for function body 3167 3168 3169@structured_monoid(tensor.arcsin) 3170def arcsin(x): 3171 """ 3172 Elemwise arcsinus of `x`. 3173 3174 """ 3175 # see decorator for function body 3176 3177 3178@structured_monoid(tensor.arctan) 3179def arctan(x): 3180 """ 3181 Elemwise arctan of `x`. 3182 3183 """ 3184 # see decorator for function body 3185 3186 3187@structured_monoid(tensor.sinh) 3188def sinh(x): 3189 """ 3190 Elemwise sinh of `x`. 3191 3192 """ 3193 # see decorator for function body 3194 3195 3196@structured_monoid(tensor.arcsinh) 3197def arcsinh(x): 3198 """ 3199 Elemwise arcsinh of `x`. 3200 3201 """ 3202 # see decorator for function body 3203 3204 3205@structured_monoid(tensor.tanh) 3206def tanh(x): 3207 """ 3208 Elemwise tanh of `x`. 3209 3210 """ 3211 # see decorator for function body 3212 3213 3214@structured_monoid(tensor.arctanh) 3215def arctanh(x): 3216 """ 3217 Elemwise arctanh of `x`. 3218 3219 """ 3220 # see decorator for function body 3221 3222 3223@structured_monoid(tensor.round_half_to_even) 3224def rint(x): 3225 """ 3226 Elemwise round half to even of `x`. 3227 3228 """ 3229 # see decorator for function body 3230 3231# Give it a simple name instead of the complex one that would automatically 3232# be derived from `tensor.round_half_to_even`. 3233rint.__name__ = 'rint' 3234 3235 3236@structured_monoid(tensor.sgn) 3237def sgn(x): 3238 """ 3239 Elemwise signe of `x`. 3240 3241 """ 3242 # see decorator for function body 3243 3244 3245@structured_monoid(tensor.ceil) 3246def ceil(x): 3247 """ 3248 Elemwise ceiling of `x`. 3249 3250 """ 3251 # see decorator for function body 3252 3253 3254@structured_monoid(tensor.floor) 3255def floor(x): 3256 """ 3257 Elemwise floor of `x`. 3258 3259 """ 3260 # see decorator for function body 3261 3262 3263@structured_monoid(tensor.log1p) 3264def log1p(x): 3265 """ 3266 Elemwise log(1 + `x`). 3267 3268 """ 3269 # see decorator for function body 3270 3271 3272@structured_monoid(tensor.expm1) 3273def expm1(x): 3274 """ 3275 Elemwise e^`x` - 1. 3276 3277 """ 3278 # see decorator for function body 3279 3280 3281@structured_monoid(tensor.deg2rad) 3282def deg2rad(x): 3283 """ 3284 Elemwise degree to radian. 3285 3286 """ 3287 # see decorator for function body 3288 3289 3290@structured_monoid(tensor.rad2deg) 3291def rad2deg(x): 3292 """ 3293 Elemwise radian to degree. 3294 3295 """ 3296 # see decorator for function body 3297 3298 3299@structured_monoid(tensor.trunc) 3300def trunc(x): 3301 """ 3302 Elemwise truncature. 3303 3304 """ 3305 # see decorator for function body 3306 3307 3308@structured_monoid(tensor.sqr) 3309def sqr(x): 3310 """ 3311 Elemwise `x` * `x`. 3312 3313 """ 3314 # see decorator for function body 3315 3316 3317@structured_monoid(tensor.sqrt) 3318def sqrt(x): 3319 """ 3320 Elemwise square root of `x`. 3321 3322 """ 3323 # see decorator for function body 3324 3325 3326@structured_monoid(tensor.conj) 3327def conj(x): 3328 """ 3329 Elemwise complex conjugate of `x`. 3330 3331 """ 3332 # see decorator for function body 3333 3334 3335class TrueDot(gof.op.Op): 3336 3337 # TODO 3338 # Simplify code by splitting into DotSS and DotSD. 3339 3340 __props__ = () 3341 # The grad_preserves_dense attribute doesn't change the 3342 # execution behavior. To let the optimizer merge nodes with 3343 # different values of this attribute we shouldn't compare it 3344 # here. 3345 3346 def __init__(self, grad_preserves_dense=True): 3347 self.grad_preserves_dense = grad_preserves_dense 3348 3349 def make_node(self, x, y): 3350 # NOTE 3351 # Because of trickiness of implementing, 3352 # we assume that the left argument x is a 3353 # SparseVariable (not dense) 3354 3355 if x.type.dtype != y.type.dtype: 3356 raise NotImplementedError() 3357 3358 if not _is_sparse_variable(x): 3359 raise TypeError(x) 3360 3361 # These are the conversions performed by scipy.sparse.dot 3362 if x.type.format == "csc" or x.type.format == "coo": 3363 myformat = "csc" 3364 elif x.type.format == "csr": 3365 myformat = "csr" 3366 else: 3367 raise NotImplementedError() 3368 3369 inputs = [x, y] # Need to convert? e.g. assparse 3370 outputs = [SparseType(dtype=x.type.dtype, format=myformat)()] 3371 return gof.Apply(self, inputs, outputs) 3372 3373 def perform(self, node, inp, out_): 3374 # TODO 3375 # -Verify that output is sufficiently sparse, 3376 # and raise a warning if it is not. 3377 # -Also determine that we are storing the 3378 # output in the best storage format? 3379 3380 x, y = inp 3381 out, = out_ 3382 rval = x.dot(y) 3383 if not scipy.sparse.issparse(rval): 3384 rval = getattr(scipy.sparse, x.format + '_matrix')(rval) 3385 # x.dot call tocsr() that will "upcast" to ['int8', 'uint8', 'short', 3386 # 'ushort', 'intc', 'uintc', 'longlong', 'ulonglong', 'single', 3387 # 'double', 'longdouble', 'csingle', 'cdouble', 'clongdouble'] 3388 # But ulonglong is uint64 on x86-64, but with a different typenum! 3389 if rval.dtype.num != np.dtype(str(rval.dtype)).num: 3390 assert str(rval.dtype) == node.outputs[0].dtype 3391 # Create a view with the expected typenum. 3392 format = node.outputs[0].type.format 3393 data = rval.data.view(dtype=node.outputs[0].dtype) 3394 indices = rval.indices 3395 indptr = rval.indptr 3396 shape = rval.shape 3397 # No need to copy indices and indptr as in CSM.perform(), 3398 # as there is only one user of them. 3399 if format == 'csc': 3400 rval = scipy.sparse.csc_matrix((data, indices, indptr), 3401 shape, copy=False) 3402 else: 3403 assert format == 'csr' 3404 rval = scipy.sparse.csr_matrix((data, indices, indptr), 3405 shape, copy=False) 3406 out[0] = rval 3407 3408 def grad(self, inputs, gout): 3409 (x, y) = inputs 3410 (gz,) = gout 3411 assert _is_sparse_variable(gz) 3412 assert _is_sparse_variable(x) 3413 3414 rval = [true_dot(gz, y.T), true_dot(x.T, gz)] 3415 if _is_dense_variable(y): 3416 if self.grad_preserves_dense: 3417 rval[1] = dense_from_sparse(rval[1]) 3418 return rval 3419 3420 def infer_shape(self, node, shapes): 3421 return [(shapes[0][0], shapes[1][1])] 3422 3423 3424def true_dot(x, y, grad_preserves_dense=True): 3425 """ 3426 Operation for efficiently calculating the dot product when 3427 one or all operands are sparse. Supported formats are CSC and CSR. 3428 The output of the operation is sparse. 3429 3430 Parameters 3431 ---------- 3432 x 3433 Sparse matrix. 3434 y 3435 Sparse matrix or 2d tensor variable. 3436 grad_preserves_dense : bool 3437 If True (default), makes the grad of dense inputs dense. 3438 Otherwise the grad is always sparse. 3439 3440 Returns 3441 ------- 3442 The dot product `x`.`y` in a sparse format. 3443 3444 Notex 3445 ----- 3446 The grad implemented is regular, i.e. not structured. 3447 3448 """ 3449 # TODO 3450 # Maybe the triple-transposition formulation 3451 # (when x is dense) is slow. See if there is a 3452 # direct way to do this. 3453 3454 if hasattr(x, 'getnnz'): 3455 x = as_sparse_variable(x) 3456 assert x.format in ["csr", "csc"] 3457 if hasattr(y, 'getnnz'): 3458 y = as_sparse_variable(y) 3459 assert y.format in ["csr", "csc"] 3460 3461 x_is_sparse_variable = _is_sparse_variable(x) 3462 y_is_sparse_variable = _is_sparse_variable(y) 3463 3464 if not x_is_sparse_variable and not y_is_sparse_variable: 3465 raise TypeError() 3466 if x_is_sparse_variable: 3467 return TrueDot(grad_preserves_dense)(x, y) 3468 else: 3469 assert y_is_sparse_variable 3470 return transpose(TrueDot(grad_preserves_dense)(y.T, x.T)) 3471 3472 3473# Dot 3474class StructuredDot(gof.Op): 3475 # See doc in instance of this Op or function after this class definition. 3476 __props__ = () 3477 3478 def make_node(self, a, b): 3479 3480 a = as_sparse_variable(a) 3481 assert a.format in ["csr", "csc", "bsr"] 3482 3483 if not _is_sparse_variable(a): 3484 raise TypeError('First argument must be of type SparseVariable ' 3485 'or SparseConstant') 3486 dtype_out = scalar.upcast(a.type.dtype, b.type.dtype) 3487 if b.type.ndim != 2: 3488 raise NotImplementedError('non-matrix b') 3489 3490 if _is_sparse_variable(b): 3491 return gof.Apply(self, [a, b], 3492 [SparseType(a.type.format, dtype_out)()]) 3493 else: 3494 return gof.Apply(self, [a, b], 3495 [tensor.tensor(dtype_out, 3496 (False, b.type.broadcastable[1]))]) 3497 3498 def perform(self, node, inputs, outputs): 3499 (a, b) = inputs 3500 (out,) = outputs 3501 if a.shape[1] != b.shape[0]: 3502 raise ValueError('shape mismatch in StructuredDot.perform', 3503 (a.shape, b.shape)) 3504 3505 variable = a * b 3506 if isinstance(node.outputs[0].type, SparseType): 3507 assert _is_sparse(variable) 3508 out[0] = variable 3509 return 3510 3511 assert _is_dense(variable) # scipy 0.7 automatically converts to dense 3512 3513 # dot of an NxM sparse matrix, with a Mx1 dense matrix, returns vector 3514 # not matrix 3515 if variable.ndim == 1: 3516 variable = np.expand_dims(variable, 1) 3517 elif variable.ndim != 2: 3518 raise Exception('Output of structured dot should be a matrix ' 3519 '(ndim=2)') 3520 3521 assert variable.ndim == 2 3522 3523 if variable.shape != (a.shape[0], b.shape[1]): 3524 if b.shape[0] == 1: 3525 raise Exception("a.shape=%s, b.shape=%s, " 3526 "variable.shape=%s ??? This is probably " 3527 "because scipy.csc_matrix dot has a bug " 3528 "with singleton dimensions (i.e. " 3529 "b.shape[0]=1), for scipy 0.6. Use scipy " 3530 "0.7. NB you have scipy version %s" % 3531 (a.shape, b.shape, variable.shape, 3532 scipy.__version__)) 3533 else: 3534 raise Exception("a.shape=%s, b.shape=%s, variable.shape=%s " 3535 " ??? I have no idea why") 3536 3537 # The cast is needed as otherwise we hit the bug mentioned into 3538 # theano._asarray function documentation. 3539 out[0] = theano._asarray(variable, str(variable.dtype)) 3540 3541 def grad(self, inputs, gout): 3542 # a is sparse, b is dense, g_out is dense 3543 # ga = g_out x b.T 3544 # gb = a.T x g_out 3545 (a, b) = inputs 3546 (g_out,) = gout 3547 return [structured_dot_grad(a, b, g_out), structured_dot(a.T, g_out)] 3548 3549 def infer_shape(self, node, shapes): 3550 return [(shapes[0][0], shapes[1][1])] 3551 3552_structured_dot = StructuredDot() 3553 3554 3555def structured_dot(x, y): 3556 """ 3557 Structured Dot is like dot, except that only the 3558 gradient wrt non-zero elements of the sparse matrix 3559 `a` are calculated and propagated. 3560 3561 The output is presumed to be a dense matrix, and is represented by a 3562 TensorType instance. 3563 3564 Parameters 3565 ---------- 3566 a 3567 A sparse matrix. 3568 b 3569 A sparse or dense matrix. 3570 3571 Returns 3572 ------- 3573 A sparse matrix 3574 The dot product of `a` and `b`. 3575 3576 Notes 3577 ----- 3578 The grad implemented is structured. 3579 3580 """ 3581 3582 # @todo: Maybe the triple-transposition formulation (when x is dense) 3583 # is slow. See if there is a direct way to do this. 3584 # (JB 20090528: Transposing tensors and sparse matrices is constant-time, 3585 # inplace, and fast.) 3586 3587 if hasattr(x, 'getnnz'): 3588 x = as_sparse_variable(x) 3589 assert x.format in ["csr", "csc"] 3590 if hasattr(y, 'getnnz'): 3591 y = as_sparse_variable(y) 3592 assert y.format in ["csr", "csc"] 3593 3594 x_is_sparse_variable = _is_sparse_variable(x) 3595 y_is_sparse_variable = _is_sparse_variable(y) 3596 if not x_is_sparse_variable and not y_is_sparse_variable: 3597 raise TypeError('structured_dot requires at least one sparse argument') 3598 3599 if x_is_sparse_variable: 3600 return _structured_dot(x, y) 3601 else: 3602 assert y_is_sparse_variable 3603 return _structured_dot(y.T, x.T).T 3604 3605 3606class StructuredDotGradCSC(gof.Op): 3607 # Op that produces the grad of StructuredDot. 3608 3609 # :param a_indices: Matrix indices 3610 # :param a_indptr: Matrix indptr 3611 # :param b: Right operand 3612 # :param g_ab: Accumulated gradient. 3613 3614 # :return: The grad of `a`.`b` for `a` accumulated 3615 # with g_ab. 3616 3617 # :note: The grad implemented is structured. 3618 # :note: a_* are the corresponding properties of a sparse 3619 # matrix in csc format. 3620 __props__ = () 3621 3622 def make_node(self, a_indices, a_indptr, b, g_ab): 3623 return gof.Apply(self, [a_indices, a_indptr, b, g_ab], 3624 [tensor.tensor(g_ab.dtype, (False,))]) 3625 3626 def perform(self, node, inputs, outputs): 3627 (a_indices, a_indptr, b, g_ab) = inputs 3628 (out,) = outputs 3629 g_a_data = np.zeros(a_indices.shape, dtype=g_ab.dtype) 3630 for j in xrange(len(a_indptr) - 1): 3631 ind0 = a_indptr[j] 3632 ind1 = a_indptr[j + 1] 3633 for i_idx in xrange(ind0, ind1): 3634 i = a_indices[i_idx] 3635 # Depending on the type of g_ab and b (sparse or dense), 3636 # the following dot product can result in a scalar or 3637 # a (1, 1) sparse matrix. 3638 dot_val = np.dot(g_ab[i], b[j].T) 3639 if isinstance(dot_val, scipy.sparse.spmatrix): 3640 dot_val = dot_val[0, 0] 3641 g_a_data[i_idx] = dot_val 3642 out[0] = g_a_data 3643 3644 def c_code_cache_version(self): 3645 return (1,) 3646 3647 def c_code(self, node, name, inputs, outputs, sub): 3648 3649 (_indices, _indptr, _d, _g) = inputs 3650 (_zout,) = outputs 3651 if node.inputs[2].type.dtype in ('complex64', 'complex128'): 3652 raise NotImplementedError('Complex types are not supported for b') 3653 if node.inputs[3].type.dtype in ('complex64', 'complex128'): 3654 raise NotImplementedError('Complex types are not supported for ' 3655 'g_ab') 3656 3657 return """ 3658 if (PyArray_NDIM(%(_d)s) != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(d) != 2"); %(fail)s;} 3659 if (PyArray_NDIM(%(_g)s) != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(g) != 2"); %(fail)s;} 3660 if (PyArray_NDIM(%(_indices)s) != 1) {PyErr_SetString(PyExc_NotImplementedError, "rank(indices) != 1"); %(fail)s;} 3661 if (PyArray_NDIM(%(_indptr)s) != 1) {PyErr_SetString(PyExc_NotImplementedError, "rank(indptr) != 1"); %(fail)s;} 3662 3663 if( PyArray_TYPE(%(_indices)s) != NPY_INT32) { 3664 PyErr_SetString(PyExc_NotImplementedError, "C"); %(fail)s;} 3665 3666 if( PyArray_TYPE(%(_indptr)s) != NPY_INT32) 3667 {PyErr_SetString(PyExc_NotImplementedError, "D"); %(fail)s;} 3668 3669 if( PyArray_DIMS(%(_d)s)[1] != PyArray_DIMS(%(_g)s)[1]) 3670 {PyErr_SetString(PyExc_NotImplementedError, "d and g have different numbers of columns"); %(fail)s;} 3671 3672 if (!%(_zout)s 3673 || (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0])) 3674 { 3675 Py_XDECREF(%(_zout)s); 3676 %(_zout)s = (PyArrayObject*) PyArray_SimpleNew(1, PyArray_DIMS(%(_indices)s), PyArray_TYPE(%(_g)s)); 3677 } 3678 3679 { //makes it compile even though labels jump over variable definitions. 3680 npy_intp nnz = PyArray_DIMS(%(_indices)s)[0]; 3681 npy_intp N = PyArray_DIMS(%(_indptr)s)[0]-1; //TODO: error checking with this 3682 3683 npy_intp Sindices = PyArray_STRIDES(%(_indices)s)[0]/PyArray_DESCR(%(_indices)s)->elsize; 3684 npy_intp Sindptr = PyArray_STRIDES(%(_indptr)s)[0]/PyArray_DESCR(%(_indptr)s)->elsize; 3685 3686 const npy_intp Sd1 = PyArray_STRIDES(%(_d)s)[1]/PyArray_DESCR(%(_d)s)->elsize; 3687 const npy_intp Sg1 = PyArray_STRIDES(%(_g)s)[1]/PyArray_DESCR(%(_g)s)->elsize; 3688 3689 const npy_intp K = PyArray_DIMS(%(_d)s)[1]; 3690 3691 const npy_int32 * __restrict__ indptr = (npy_int32 *)PyArray_DATA(%(_indptr)s); 3692 const npy_int32 * __restrict__ indices = (npy_int32 *)PyArray_DATA(%(_indices)s); 3693 3694 // loop over columns 3695 for (npy_int32 j = 0; j < N; ++j) 3696 { 3697 // extract j-th row of dense matrix 3698 const dtype_%(_d)s* __restrict__ d_row = (dtype_%(_d)s*)(PyArray_BYTES(%(_d)s) + PyArray_STRIDES(%(_d)s)[0] * j); 3699 if(j >= PyArray_DIMS(%(_d)s)[0]) {PyErr_SetString(PyExc_NotImplementedError, "G"); %(fail)s;} 3700 3701 // for each non-null value in the sparse column 3702 for (npy_int32 i_idx = indptr[j * Sindptr]; i_idx < indptr[(j+1) * Sindptr]; ++i_idx) 3703 { 3704 // extract row index of non-null value 3705 npy_int32 i = indices[i_idx * Sindices]; 3706 3707 // extract corresponding row in gradient 3708 const dtype_%(_g)s* __restrict__ g_row = (dtype_%(_g)s*)(PyArray_BYTES(%(_g)s) + PyArray_STRIDES(%(_g)s)[0] * i); 3709 double ip = 0.0; 3710 3711 // make sure that row index is not bigger than actual number of rows 3712 // Note: wouldn't the above operation fail if that were the case ? 3713 // when would this ever be true anyway ? 3714 if (i >= PyArray_DIMS(%(_g)s)[0]) 3715 {PyErr_SetString(PyExc_NotImplementedError, "H"); %(fail)s;} 3716 3717 // perform dot product of dense and sparse rows 3718 for(int k = 0; k < K; ++k) 3719 { 3720 ip += d_row[k * Sd1] * g_row[k*Sg1]; 3721 } 3722 3723 // write resulting gradient to sparse output 3724 ((dtype_%(_zout)s* __restrict__)(PyArray_BYTES(%(_zout)s) + i_idx * PyArray_STRIDES(%(_zout)s)[0]))[0] = ip; 3725 } 3726 } 3727 } 3728 3729 """ % dict(locals(), **sub) 3730 3731 def infer_shape(self, node, shapes): 3732 return [shapes[0]] 3733sdg_csc = StructuredDotGradCSC() 3734 3735 3736class StructuredDotGradCSR(gof.Op): 3737 # Op that produces the grad of StructuredDot. 3738 3739 # :param a_indices: Matrix indices 3740 # :param a_indptr: Matrix indptr 3741 # :param b: Right operand 3742 # :param g_ab: Accumulated gradient. 3743 3744 # :return: The grad of `a`.`b` for `a` accumulated 3745 # with g_ab. 3746 3747 # :note: The grad implemented is structured. 3748 # :note: a_* are the corresponding properties of a sparse 3749 # matrix in csr format. 3750 __props__ = () 3751 3752 def make_node(self, a_indices, a_indptr, b, g_ab): 3753 return gof.Apply(self, [a_indices, a_indptr, b, g_ab], 3754 [tensor.tensor(b.dtype, (False,))]) 3755 3756 def perform(self, node, inputs, outputs): 3757 (a_indices, a_indptr, b, g_ab) = inputs 3758 (out,) = outputs 3759 g_a_data = np.zeros(a_indices.shape, dtype=g_ab.dtype) 3760 for i in xrange(len(a_indptr) - 1): # loop over rows 3761 ind0 = a_indptr[i] 3762 ind1 = a_indptr[i + 1] 3763 # loop over values in that row (columns) 3764 for j_idx in xrange(ind0, ind1): 3765 j = a_indices[j_idx] 3766 # grad is dot product of i-th row of gradient with j-th row of b 3767 # Depending on the type of g_ab and b (sparse or dense), 3768 # the following dot product can result in a scalar or 3769 # a (1, 1) sparse matrix. 3770 dot_val = np.dot(g_ab[i], b[j].T) 3771 if isinstance(dot_val, scipy.sparse.spmatrix): 3772 dot_val = dot_val[0, 0] 3773 g_a_data[j_idx] = dot_val 3774 out[0] = g_a_data 3775 3776 def c_code_cache_version(self): 3777 return (1,) 3778 3779 def c_code(self, node, name, inputs, outputs, sub): 3780 3781 (_indices, _indptr, _d, _g) = inputs 3782 (_zout,) = outputs 3783 if node.inputs[2].type.dtype in ('complex64', 'complex128'): 3784 raise NotImplementedError('Complex types are not supported for b') 3785 if node.inputs[3].type.dtype in ('complex64', 'complex128'): 3786 raise NotImplementedError('Complex types are not supported for ' 3787 'g_ab') 3788 3789 return """ 3790 if (PyArray_NDIM(%(_d)s) != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(d) != 2"); %(fail)s;} 3791 if (PyArray_NDIM(%(_g)s) != 2) {PyErr_SetString(PyExc_NotImplementedError, "rank(g) != 2"); %(fail)s;} 3792 if (PyArray_NDIM(%(_indices)s) != 1) {PyErr_SetString(PyExc_NotImplementedError, "rank(indices) != 1"); %(fail)s;} 3793 if (PyArray_NDIM(%(_indptr)s) != 1) {PyErr_SetString(PyExc_NotImplementedError, "rank(indptr) != 1"); %(fail)s;} 3794 3795 if( PyArray_TYPE(%(_indices)s) != NPY_INT32) { 3796 PyErr_SetString(PyExc_NotImplementedError, "C"); %(fail)s;} 3797 3798 if( PyArray_TYPE(%(_indptr)s) != NPY_INT32) 3799 {PyErr_SetString(PyExc_NotImplementedError, "D"); %(fail)s;} 3800 3801 if( PyArray_DIMS(%(_d)s)[1] != PyArray_DIMS(%(_g)s)[1]) 3802 {PyErr_SetString(PyExc_NotImplementedError, "d and g have different numbers of columns"); %(fail)s;} 3803 3804 if (!%(_zout)s 3805 || (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0])) 3806 { 3807 Py_XDECREF(%(_zout)s); 3808 %(_zout)s = (PyArrayObject*) PyArray_SimpleNew(1, PyArray_DIMS(%(_indices)s), PyArray_TYPE(%(_g)s)); 3809 } 3810 3811 { //makes it compile even though labels jump over variable definitions. 3812 npy_intp nnz = PyArray_DIMS(%(_indices)s)[0]; 3813 // extract number of rows 3814 npy_intp N = PyArray_DIMS(%(_indptr)s)[0]-1; //TODO: error checking with this 3815 3816 npy_intp Sindices = PyArray_STRIDES(%(_indices)s)[0]/PyArray_DESCR(%(_indices)s)->elsize; 3817 npy_intp Sindptr = PyArray_STRIDES(%(_indptr)s)[0]/PyArray_DESCR(%(_indptr)s)->elsize; 3818 3819 const npy_intp Sd1 = PyArray_STRIDES(%(_d)s)[1]/PyArray_DESCR(%(_d)s)->elsize; 3820 const npy_intp Sg1 = PyArray_STRIDES(%(_g)s)[1]/PyArray_DESCR(%(_g)s)->elsize; 3821 3822 const npy_intp K = PyArray_DIMS(%(_d)s)[1]; 3823 3824 const npy_int32 * __restrict__ indptr = (npy_int32 *)PyArray_DATA(%(_indptr)s); 3825 const npy_int32 * __restrict__ indices = (npy_int32 *)PyArray_DATA(%(_indices)s); 3826 3827 // loop over columns of sparse matrix 3828 for (npy_int32 i = 0; i < N; ++i) 3829 { 3830 // for each non-null value in the sparse row 3831 for (npy_int32 j_idx = indptr[i * Sindptr]; j_idx < indptr[(i+1) * Sindptr]; ++j_idx) 3832 { 3833 // extract column index of non-null value 3834 npy_int32 j = indices[j_idx * Sindices]; 3835 3836 // extract j-th row of dense matrix 3837 const dtype_%(_d)s* __restrict__ d_row = (dtype_%(_d)s*)(PyArray_BYTES(%(_d)s) + PyArray_STRIDES(%(_d)s)[0] * j); 3838 if(j >= PyArray_DIMS(%(_d)s)[0]) {PyErr_SetString(PyExc_NotImplementedError, "G"); %(fail)s;} 3839 3840 // extract corresponding row in gradient 3841 const dtype_%(_g)s* __restrict__ g_row = (dtype_%(_g)s*)(PyArray_BYTES(%(_g)s) + PyArray_STRIDES(%(_g)s)[0] * i); 3842 double ip = 0.0; 3843 3844 // make sure that row index is not bigger than actual number of rows 3845 // Note: wouldn't the above operation fail if that were the case ? 3846 // when would this ever be true anyway ? 3847 if (i >= PyArray_DIMS(%(_g)s)[0]) 3848 {PyErr_SetString(PyExc_NotImplementedError, "H"); %(fail)s;} 3849 3850 // perform dot product of dense and sparse rows 3851 for(int k = 0; k < K; ++k) 3852 { 3853 ip += d_row[k * Sd1] * g_row[k*Sg1]; 3854 } 3855 3856 // write resulting gradient to sparse output 3857 ((dtype_%(_zout)s* __restrict__)(PyArray_BYTES(%(_zout)s) + j_idx * PyArray_STRIDES(%(_zout)s)[0]))[0] = ip; 3858 } 3859 } 3860 } 3861 3862 """ % dict(locals(), **sub) 3863 3864 def infer_shape(self, node, shapes): 3865 return [shapes[0]] 3866sdg_csr = StructuredDotGradCSR() 3867 3868 3869def structured_dot_grad(sparse_A, dense_B, ga): 3870 if sparse_A.type.format in ('csc', 'csr'): 3871 3872 if sparse_A.type.format == 'csc': 3873 sdgcsx = sdg_csc 3874 CSx = CSC 3875 else: 3876 sdgcsx = sdg_csr 3877 CSx = CSR 3878 3879 g_A_data = sdgcsx(csm_indices(sparse_A), 3880 csm_indptr(sparse_A), dense_B, ga) 3881 return CSx(g_A_data, csm_indices(sparse_A), 3882 csm_indptr(sparse_A), csm_shape(sparse_A)) 3883 else: 3884 raise NotImplementedError() 3885 3886 3887class SamplingDot(gof.op.Op): 3888 # See doc in instance of this Op or function after this class definition. 3889 __props__ = () 3890 3891 def make_node(self, x, y, p): 3892 x = tensor.as_tensor_variable(x) 3893 y = tensor.as_tensor_variable(y) 3894 p = as_sparse_variable(p) 3895 assert p.format in ["csr", "csc"] 3896 3897 if not _is_sparse_variable(p): 3898 raise TypeError(p) 3899 3900 # TODO: use it. 3901 dtype_out = scalar.upcast(x.type.dtype, y.type.dtype, p.type.dtype) # noqa 3902 3903 return gof.Apply(self, [x, y, p], [p.type()]) 3904 3905 def perform(self, node, inputs, outputs): 3906 (x, y, p) = inputs 3907 (out,) = outputs 3908 if _is_sparse(x): 3909 raise TypeError(x) 3910 3911 if _is_sparse(y): 3912 raise TypeError(y) 3913 3914 if not _is_sparse(p): 3915 raise TypeError(p) 3916 3917 out[0] = p.__class__(p.multiply(np.dot(x, y.T))) 3918 3919 def grad(self, inputs, gout): 3920 (x, y, p) = inputs 3921 (gz,) = gout 3922 rval = [ 3923 dot(p * gz, y), 3924 dot((p * gz).T, x), 3925 grad_not_implemented(self, 2, p) 3926 ] 3927 3928 return rval 3929 3930 def infer_shape(self, node, ins_shapes): 3931 return [ins_shapes[2]] 3932 3933sampling_dot = SamplingDot() 3934""" 3935Operand for calculating the dot product dot(`x`, `y`.T) = `z` when you 3936only want to calculate a subset of `z`. 3937 3938It is equivalent to `p` o (`x` . `y`.T) where o is the element-wise 3939product, `x` and `y` operands of the dot product and `p` is a matrix that 3940contains 1 when the corresponding element of `z` should be calculated 3941and 0 when it shouldn't. Note that SamplingDot has a different interface 3942than `dot` because SamplingDot requires `x` to be a `m`x`k` matrix while 3943`y` is a `n`x`k` matrix instead of the usual `k`x`n` matrix. 3944 3945Notes 3946----- 3947It will work if the pattern is not binary value, but if the 3948pattern doesn't have a high sparsity proportion it will be slower 3949then a more optimized dot followed by a normal elemwise 3950multiplication. 3951 3952The grad implemented is regular, i.e. not structured. 3953 3954Parameters 3955---------- 3956x 3957 Tensor matrix. 3958y 3959 Tensor matrix. 3960p 3961 Sparse matrix in csr format. 3962 3963Returns 3964------- 3965sparse matrix 3966 A dense matrix containing the dot product of `x` by `y`.T only 3967 where `p` is 1. 3968 3969""" 3970 3971 3972class Dot(gof.op.Op): 3973 # See doc in instance of this Op or function after this class definition. 3974 __props__ = () 3975 3976 def __str__(self): 3977 return "Sparse" + self.__class__.__name__ 3978 3979 def infer_shape(self, node, shapes): 3980 xshp, yshp = shapes 3981 x, y = node.inputs 3982 if x.ndim == 2 and y.ndim == 2: 3983 return [(xshp[0], yshp[1])] 3984 if x.ndim == 1 and y.ndim == 2: 3985 return [(yshp[1],)] 3986 if x.ndim == 2 and y.ndim == 1: 3987 return [(xshp[0],)] 3988 if x.ndim == 1 and y.ndim == 1: 3989 return [()] 3990 raise NotImplementedError() 3991 3992 def make_node(self, x, y): 3993 dtype_out = scalar.upcast(x.dtype, y.dtype) 3994 3995 # Sparse dot product should have at least one sparse variable 3996 # as input. If the other one is not sparse, it has to be converted 3997 # into a tensor. 3998 if isinstance(x, scipy.sparse.spmatrix): 3999 x = as_sparse_variable(x) 4000 if isinstance(y, scipy.sparse.spmatrix): 4001 y = as_sparse_variable(y) 4002 x_is_sparse_var = _is_sparse_variable(x) 4003 y_is_sparse_var = _is_sparse_variable(y) 4004 4005 if not x_is_sparse_var and not y_is_sparse_var: 4006 raise TypeError( 4007 "Sparse dot product should have at least one " 4008 "sparse variable as inputs, but the inputs are " 4009 "%s (%s) and %s (%s)." % (x, x.type, y, y.type)) 4010 4011 if x_is_sparse_var: 4012 broadcast_x = (False,) * x.ndim 4013 else: 4014 x = tensor.as_tensor_variable(x) 4015 broadcast_x = x.type.broadcastable 4016 assert y.format in ["csr", "csc"] 4017 if x.ndim not in (1, 2): 4018 raise TypeError( 4019 'theano.sparse.Dot: input 0 (0-indexed) must have ndim of ' 4020 '1 or 2, %d given.' % x.ndim) 4021 4022 if y_is_sparse_var: 4023 broadcast_y = (False,) * y.ndim 4024 else: 4025 y = tensor.as_tensor_variable(y) 4026 broadcast_y = y.type.broadcastable 4027 assert x.format in ["csr", "csc"] 4028 if y.ndim not in (1, 2): 4029 raise TypeError( 4030 'theano.sparse.Dot: input 1 (1-indexed) must have ndim of ' 4031 '1 or 2, %d given.' % y.ndim) 4032 4033 if len(broadcast_y) == 2: 4034 broadcast_out = broadcast_x[:-1] + broadcast_y[1:] 4035 elif len(broadcast_y) == 1: 4036 broadcast_out = broadcast_x[:-1] 4037 return gof.Apply(self, [x, y], [tensor.tensor(dtype=dtype_out, 4038 broadcastable=broadcast_out)]) 4039 4040 def perform(self, node, inputs, out): 4041 x, y = inputs 4042 out = out[0] 4043 x_is_sparse = _is_sparse(x) 4044 y_is_sparse = _is_sparse(y) 4045 4046 if not x_is_sparse and not y_is_sparse: 4047 raise TypeError(x) 4048 4049 rval = x * y 4050 4051 if x_is_sparse and y_is_sparse: 4052 rval = rval.toarray() 4053 4054 out[0] = theano._asarray(rval, dtype=node.outputs[0].dtype) 4055 4056 def grad(self, inputs, gout): 4057 (x, y) = inputs 4058 (gz,) = gout 4059 assert _is_sparse_variable(x) or _is_sparse_variable(y) 4060 rval = [] 4061 4062 if _is_dense_variable(y): 4063 rval.append(tensor.dot(gz, y.T)) 4064 else: 4065 rval.append(dot(gz, y.T)) 4066 if _is_dense_variable(x): 4067 rval.append(tensor.dot(x.T, gz)) 4068 else: 4069 rval.append(dot(x.T, gz)) 4070 4071 return rval 4072_dot = Dot() 4073 4074 4075def dot(x, y): 4076 """ 4077 Operation for efficiently calculating the dot product when 4078 one or all operands is sparse. Supported format are CSC and CSR. 4079 The output of the operation is dense. 4080 4081 Parameters 4082 ---------- 4083 x 4084 Sparse or dense matrix variable. 4085 y 4086 Sparse or dense matrix variable. 4087 4088 Returns 4089 ------- 4090 The dot product `x`.`y` in a dense format. 4091 4092 Notes 4093 ----- 4094 The grad implemented is regular, i.e. not structured. 4095 4096 At least one of `x` or `y` must be a sparse matrix. 4097 4098 When the operation has the form dot(csr_matrix, dense) 4099 the gradient of this operation can be performed inplace 4100 by UsmmCscDense. This leads to significant speed-ups. 4101 4102 """ 4103 4104 if hasattr(x, 'getnnz'): 4105 x = as_sparse_variable(x) 4106 if hasattr(y, 'getnnz'): 4107 y = as_sparse_variable(y) 4108 4109 x_is_sparse_variable = _is_sparse_variable(x) 4110 y_is_sparse_variable = _is_sparse_variable(y) 4111 4112 if not x_is_sparse_variable and not y_is_sparse_variable: 4113 raise TypeError() 4114 4115 return _dot(x, y) 4116 4117 4118class Usmm(gof.op.Op): 4119 # See doc in instance of this Op or function after this class definition. 4120 # We don't implement the infer_shape as it is 4121 # inserted by optimization only. 4122 __props__ = () 4123 4124 def __str__(self): 4125 return 'Usmm{no_inplace}' 4126 4127 def make_node(self, alpha, x, y, z): 4128 if not _is_sparse_variable(x) and not _is_sparse_variable(y): 4129 # If x and y are tensor, we don't want to use this class 4130 # We should use Dot22 and Gemm in that case. 4131 raise TypeError(x) 4132 4133 dtype_out = scalar.upcast(alpha.type.dtype, x.type.dtype, 4134 y.type.dtype, z.type.dtype) 4135 alpha = tensor.as_tensor_variable(alpha) 4136 z = tensor.as_tensor_variable(z) 4137 4138 assert z.ndim == 2 4139 assert alpha.type.broadcastable == (True,) * alpha.ndim 4140 if not _is_sparse_variable(x): 4141 x = tensor.as_tensor_variable(x) 4142 assert y.format in ["csr", "csc"] 4143 assert x.ndim == 2 4144 if not _is_sparse_variable(y): 4145 y = tensor.as_tensor_variable(y) 4146 assert x.format in ["csr", "csc"] 4147 assert y.ndim == 2 4148 4149 return gof.Apply(self, [alpha, x, y, z], 4150 [tensor.tensor(dtype=dtype_out, 4151 broadcastable=(False, False))]) 4152 4153 def perform(self, node, inputs, outputs): 4154 (alpha, x, y, z) = inputs 4155 (out,) = outputs 4156 x_is_sparse = _is_sparse(x) 4157 y_is_sparse = _is_sparse(y) 4158 4159 if not x_is_sparse and not y_is_sparse: 4160 raise TypeError(x) 4161 4162 rval = x * y 4163 if isinstance(rval, scipy.sparse.spmatrix): 4164 rval = rval.toarray() 4165 if rval.dtype == alpha.dtype: 4166 rval *= alpha # Faster because operation is inplace 4167 else: 4168 rval = rval * alpha 4169 if rval.dtype == z.dtype: 4170 rval += z # Faster because operation is inplace 4171 else: 4172 rval = rval + z 4173 4174 out[0] = rval 4175usmm = Usmm() 4176""" 4177Performs the expression `alpha` * `x` `y` + `z`. 4178 4179Parameters 4180---------- 4181x 4182 Matrix variable. 4183y 4184 Matrix variable. 4185z 4186 Dense matrix. 4187alpha 4188 A tensor scalar. 4189 4190Returns 4191------- 4192The dense matrix resulting from `alpha` * `x` `y` + `z`. 4193 4194Notes 4195----- 4196The grad is not implemented for this op. 4197At least one of `x` or `y` must be a sparse matrix. 4198 4199""" 4200 4201 4202class ConstructSparseFromList(gof.Op): 4203 # See doc in instance of this Op or function after this class definition. 4204 __props__ = () 4205 4206 def make_node(self, x, values, ilist): 4207 """ 4208 4209 Parameters 4210 ---------- 4211 x 4212 A dense matrix that specify the output shape. 4213 values 4214 A dense matrix with the values to use for output. 4215 ilist 4216 A dense vector with the same length as the number of rows of values. 4217 It specify where in the output to put the corresponding rows. 4218 4219 This create a sparse matrix with the same shape as `x`. Its 4220 values are the rows of `values` moved. Pseudo-code:: 4221 4222 output = csc_matrix.zeros_like(x, dtype=values.dtype) 4223 for in_idx, out_idx in enumerate(ilist): 4224 output[out_idx] = values[in_idx] 4225 4226 """ 4227 x_ = theano.tensor.as_tensor_variable(x) 4228 values_ = theano.tensor.as_tensor_variable(values) 4229 ilist_ = theano.tensor.as_tensor_variable(ilist) 4230 4231 if ilist_.type.dtype not in integer_dtypes: 4232 raise TypeError('index must be integers') 4233 if ilist_.type.ndim != 1: 4234 raise TypeError('index must be vector') 4235 if x_.type.ndim != 2: 4236 raise TypeError( 4237 'cannot create a sparse matrix with %d dimensions' % 4238 x_.type.ndim) 4239 if values_.type.ndim != 2: 4240 raise TypeError( 4241 'cannot create a sparse matrix from values with %d ndim' % 4242 values_.type.ndim) 4243 4244 # We only need the shape of `x` in the perform 4245 # If we keep in the graph the x variable as input of the Apply node, 4246 # this can rise the memory usage. That is why the Apply node 4247 # take `x_.shape` as input and not `x`. 4248 return gof.Apply(self, [x_.shape, values_, ilist_], 4249 [csc_matrix(dtype=x.dtype)]) 4250 4251 def perform(self, node, inp, out_): 4252 out_shape, values, ilist = inp 4253 out, = out_ 4254 rows, cols = values.shape 4255 assert rows == len(ilist) 4256 indptr = np.arange(cols + 1) * rows 4257 indices = as_strided(ilist, 4258 strides=(0, ilist.strides[0]), 4259 shape=(cols, ilist.shape[0])).flatten() 4260 data = values.T.flatten() 4261 out[0] = scipy.sparse.csc_matrix((data, indices, indptr), 4262 shape=out_shape, 4263 dtype=values.dtype) 4264 4265 def infer_shape(self, node, ishapes): 4266 x = node.inputs[0] 4267 return [[x[0], x[1]]] 4268 4269 def R_op(self, inputs, eval_points): 4270 if None in eval_points[:2]: 4271 return [None] 4272 return self.make_node(eval_points[0], eval_points[1], 4273 *inputs[2:]).outputs 4274 4275 def connection_pattern(self, node): 4276 4277 rval = [[True], [True], [False]] 4278 return rval 4279 4280 def grad(self, inputs, grads): 4281 g_output, = grads 4282 x, y = inputs[:2] 4283 idx_list = inputs[2:] 4284 4285 gx = g_output 4286 gy = theano.tensor.advanced_subtensor1(g_output, *idx_list) 4287 4288 return [gx, gy] + [DisconnectedType()()] * len(idx_list) 4289 4290construct_sparse_from_list = ConstructSparseFromList() 4291""" 4292Constructs a sparse matrix out of a list of 2-D matrix rows. 4293 4294Notes 4295----- 4296The grad implemented is regular, i.e. not structured. 4297 4298""" 4299