[SciPy-User] loadmat/savemat sparse with umfpack problem

[Robin]

On Mon, Dec 14, 2009 at 8:29 AM, Robert Cimrman <cimrman3 at ntc.zcu.cz> wrote:
> Hi Robin,
>
> Robin wrote:
>> Hi,
>>
>> I have a problem with sparse matrices that have been roundtripped
>> through savemat/loadmat and scikits.umf:
>> In [29]: scipy.__version__
>> Out[29]: '0.8.0.dev6136'
>>
>> scikits.umfpack doesn't have a version but it is latest from scikits
>> svn (r2239).
>>
>> Here is the simplest way I could recreate the problem:
>>
>> from scikits.umfpack import UmfpackContext
>> import scipy.sparse as sparse
>> from scipy.io import loadmat, savemat
>>
>> a = sparse.eye(3,3,format='csc')
>>
>> umf = UmfpackContext()
>> print 'Original sparse matrix:'
>> print a.__repr__()
>> # works fine
>> umf.numeric(a)
>> print 'savemat/loadmat ...'
>> savemat('test',{'a':a})
>> a2 = loadmat('test')['a']
>> print 'Loaded sparse matrix:'
>> print a2.__repr__()
>> # doesnt work
>> umf.numeric(a2)
>>
>> which outputs:
>>
>> Original sparse matrix:
>> <3x3 sparse matrix of type '<type 'numpy.float64'>'
>>         with 3 stored elements in Compressed Sparse Column format>
>> savemat/loadmat ...
>> Loaded sparse matrix:
>> <3x3 sparse matrix of type '<type 'numpy.float64'>'
>>         with 3 stored elements in Compressed Sparse Column format>
>> ---------------------------------------------------------------------------
>> TypeError                                 Traceback (most recent call last)
>>
>> /Users/robince/svn/pyentropy/pyentropy/umf.py in <module>()
>>      15 print 'Loaded sparse matrix:'
>>      16 print a2.__repr__()
>> ---> 17 umf.numeric(a2)
>>      18
>>      19
>>
>> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/scikits.umfpack-5.1.0-py2.5-macosx-10.3-i386.egg/scikits/umfpack/umfpack.pyc
>> in numeric(self, mtx)
>>     393
>>     394         if self._symbolic is None:
>> --> 395             self.symbolic( mtx )
>>     396
>>     397         indx = self._getIndx( mtx )
>>
>> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/scikits.umfpack-5.1.0-py2.5-macosx-10.3-i386.egg/scikits/umfpack/umfpack.pyc
>> in symbolic(self, mtx)
>>     364                     = self.funs.symbolic( mtx.shape[0], mtx.shape[1],
>>     365                                           mtx.indptr, indx, mtx.data,
>> --> 366                                           self.control, self.info )
>>     367         else:
>>     368             real, imag = mtx.data.real.copy(), mtx.data.imag.copy()
>>
>> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/scikits.umfpack-5.1.0-py2.5-macosx-10.3-i386.egg/scikits/umfpack/_umfpack.pyc
>> in umfpack_di_symbolic(*args)
>>     435         double Control, double Info) -> int
>>     436     """
>> --> 437   return __umfpack.umfpack_di_symbolic(*args)
>>     438
>>     439 def umfpack_dl_symbolic(*args):
>>
>> TypeError: not a C array
>> WARNING: Failure executing file: <umf.py>
>>
>> I can't figure out whats causing it - and whether its a bug in
>> savemat/loadmat or scikits.umfpack.
>>
>> Is scikits.umfpack still supported? Is there a way to prefactor matrix
>> with the libraries built into scipy (I need to solve the same large
>> sparse matrix many times so I was prefactoring with umfpack).
>
> I am not sure what causes this, but you try the following:
>
> 1. umfpack requires the indices array to be sorted in ascending order, there is
> a function to ensure that:
>
> csc_matrix.ensure_sorted_indices()
>
> 2. the wrappers expect all the sparse matrix arrays (indptr, indices, data) to
> be in c-contiguous order - try a2.indices = a2.indices.copy() etc.
>
> Hope that helps,
> r.

Hi,

Thanks very much. Sorting the indices seems to fix it. (although
ensure_sorted_indices is deprecated for sorted_indices or
sort_indices).

Still not sure why though - in this simple example the indices do seem
to be sorted - the only change I can find from calling the
sorted_indices function is that after the WRITEABLE flag of .indices
is True (as it is originally) but on the loadmat'ed array WRITEABLE is
False.

Could this be causing it? Is it a bug in loadmat - should the flag be
set differently? (I'm not sure what it does).

In [31]: a.indices   # origina
Out[31]: array([0, 1, 2])
In [32]: a2.indices  # loaded
Out[32]: array([0, 1, 2])
In [33]: a2.sorted_indices().indices  # loaded + sorted
Out[33]: array([0, 1, 2])
In [34]: a.indices.flags # original
Out[34]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
In [35]: a2.indices.flags # laoded
Out[35]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : False
  ALIGNED : True
  UPDATEIFCOPY : False
In [37]: a2.sorted_indices().indices.flags  # loaded and sorted
Out[37]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

Thanks again for the help,

Cheers

Robin