docs: def unbroadcast(a, b):
      '''
      unbroadcast(a, b) yields a tuple (aa, bb) that is equivalent to (a, b) except that aa and bb
        have been reshaped such that arithmetic numpy operations such as aa * bb will result in
        row-wise operation instead of column-wise broadcasting.
      '''
      # they could be sparse:
      spa = sps.issparse(a)
      spb = sps.issparse(b)
      if   spa and spb: return (a,b)
      elif spa or  spb:
          def fix(sp,nm):
              nm = np.asarray(nm)
              dnm = len(nm.shape)
              nnm = np.prod(nm.shape)
              # if we have (sparse matrix) * (high-dim array), unbroadcast the dense array
              if   dnm == 0: return (sp, np.reshape(nm, (1,   1)))
              elif dnm == 1: return (sp, np.reshape(nm, (nnm, 1)))
              elif dnm == 2: return (sp, nm)
              else:          return unbroadcast(sp.toarray(), nm)
          return fix(a, b) if spa else tuple(reversed(fix(b, a)))
      # okay, no sparse matrices found:
      a = np.asarray(a)
      b = np.asarray(b)
      da = len(a.shape)
      db = len(b.shape)
      if   da > db: return (a, np.reshape(b, b.shape + tuple(np.ones(da-db, dtype=np.int))))
      elif da < db: return (np.reshape(a, a.shape + tuple(np.ones(db-da, dtype=np.int))), b)
      else:         return (a, b)