Thread Subject: spfun

Hello, I was wondering if anyone knows a solution to divide element-
wise two sparse matrices, A./B, but considering only the positions
where A is non-zero.

I though in using spfun, but it works only with one variable input.

Thank you in advance.

Hi Arl,

what about this?

C = A.*spfun(@(x) 1./x,B);

Best regards, anon

arl <arlourenco@gmail.com> wrote in message <f495641d-048f-423d-a618-676f969fa814@p36g2000vbn.googlegroups.com>...
> Hello, I was wondering if anyone knows a solution to divide element-
> wise two sparse matrices, A./B, but considering only the positions
> where A is non-zero.
>
> I though in using spfun, but it works only with one variable input.
>
> Thank you in advance.

arl <arlourenco@gmail.com> wrote in message <f495641d-048f-423d-a618-676f969fa814@p36g2000vbn.googlegroups.com>...
> Hello, I was wondering if anyone knows a solution to divide element-
> wise two sparse matrices, A./B, but considering only the positions
> where A is non-zero.
>
> I though in using spfun, but it works only with one variable input.
>
> Thank you in advance.

result=A;
idx=logical(A);
result(idx)=A(idx)./B(idx);

"Matt " <xys@whatever.com> wrote in message <heedra$cuq$1@fred.mathworks.com>...

>
> result=A;
> idx=logical(A);
> result(idx)=A(idx)./B(idx);

CAREFUL!!!! Matlab is known for being BUGGY in doing such indexing with sparse

result(idx)=A(idx)./B(idx);

>> A=sparse(100000,100000,1)

A =

             (100000,100000) 1

>> B=A

B =

             (100000,100000) 1

>> result=A;
>> idx=logical(A);
>> result(idx)=A(idx)./B(idx)

result =

              (65408,14101) 1
             (100000,100000) 1

%%%%%%%%%%%%%%%%%%

This is one of the reason why I wrote this package
http://www.mathworks.com/matlabcentral/fileexchange/23488-sparse-sub-access

>> [i j]=find(A);
>> AB = setsparse(A,i,j,getsparse(B,i,j),@(a,b) a./b)

AB =

             (100000,100000) 1

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <heefss$lcl$1@fred.mathworks.com>...

> %%%%%%%%%%%%%%%%%%
>
> This is one of the reason why I wrote this package
> http://www.mathworks.com/matlabcentral/fileexchange/23488-sparse-sub-access
>
> >> [i j]=find(A);
> >> AB = setsparse(A,i,j,getsparse(B,i,j),@(a,b) a./b)
>
> AB =
>
> (100000,100000) 1
>
> % Bruno


I didn't realize the extent of the problem. It wonder if one could provide a subclass of @double to interface to this package. That way, we could use the old indexing syntax instead of things like this

AB = setsparse(A,i,j,getsparse(B,i,j),@(a,b) a./b)

"Matt " <xys@whatever.com> wrote in message <heeinc$ha7$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <heefss$lcl$1@fred.mathworks.com>...
>
> > %%%%%%%%%%%%%%%%%%
> >
> > This is one of the reason why I wrote this package
> > http://www.mathworks.com/matlabcentral/fileexchange/23488-sparse-sub-access
> >
> > >> [i j]=find(A);
> > >> AB = setsparse(A,i,j,getsparse(B,i,j),@(a,b) a./b)
> >
> > AB =
> >
> > (100000,100000) 1
> >
> > % Bruno
>
>
> I didn't realize the extent of the problem. It wonder if one could provide a subclass of @double to interface to this package. That way, we could use the old indexing syntax instead of things like this

Could be a nice project indeed, if there is demand out there. I'm not going to involve in coding it though (lack of time).

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <heejqv$q09$1@fred.mathworks.com>...
>
> Could be a nice project indeed, if there is demand out there. I'm not going to involve in coding it though (lack of time).
>
=======================


Would you be willing to at least test-drive the class I've coded at the bottom of this post (to be put in a file called RobustSparse.m). It will take me some time to fully understand/install your package, whereas my code would be a 1-file installation for you.

If you're willing, the following revised tests should give the correct results. Otherwise, it will default with a warning to the normal sparse indexing behavior

>> A=RobustSparse(100000,100000,1), B=A
 
B =
 
             (100000,100000) 1

>> result=A;
>> idx=logical(A);
>> result(idx)=A(idx)./B(idx)





%%%%Put in file RobustSparse.m

classdef RobustSparse < double
%RobustSparse - a class essentially the same as sparse @double, but with
%some bug protection in subsref/subsasgn courtesy of Bruno Luong.

    methods
       
        function obj=RobustSparse(varargin)
            %constructor for RobustSparse

            data=sparse(varargin{:});
            obj=obj@double(data);
            
        end
        
        function objnew=subsref(obj,S)
        %SUBSREF for RobustSparse class - use Bruno Luong's routines
        
            nn=length(S.subs);
            idx=S.subs{1};
            
            if nn==1 && ~ischar(idx)
                

                try


                   if islogical(idx)
                       
                       [i,j]=find(idx);
                       
                   else%linear index
                       
                       [i,j]=ind2sub(size(obj),idx);
                       
                   end
                   
                  objnew=getsparse(double(obj),i,j);
                  objnew=reshape(objnew,size(idx));

                 


                catch%Default to usual method

                    disp 'sparse-sub-access not available. Defaulting...'
                    
                   objnew=subsref@double(obj,S);

                end

            else
                
                objnew=subsref@double(obj,S);
                
            end
            
           objnew=RobustSparse(objnew);
           
        end
        
        function objnew=subsasgn(obj,S,rhs)
        %SUBSASGN for RobustSparse class - use Bruno Luong's routines
            
            nn=length(S.subs);
            idx=S.subs{1};
            
            if nn==1 && ~ischar(idx)
                

                try


                   if islogical(idx)
                       
                       [i,j]=find(idx);
                       
                   else%single linear index
                       
                       [i,j]=ind2sub(size(obj),idx);
                       
                   end
                   
                  objnew=setsparse(double(obj),i,j,rhs);

                 


                catch%Default to usual method

                    disp 'sparse-sub-access not available. Defaulting...'
                    
                   objnew=subsasgn@double(obj,S,rhs);

                end

            else
                
                objnew=subsasgn@double(obj,S,rhs);
                
            end
        
           objnew=RobustSparse(objnew);
        
        

        end
        

        
        function obj=uplus(obj)
            
        end

        function obj=uminus(obj)
            
            obj=RobustSparse(uminus@double(obj));
            
        end
        
        function obj=plus(L,R)
            
          obj = RobustSparse( plus@double(L,R) );
            
        end

        function obj=minus(L,R)
            
          obj = RobustSparse( minus@double(L,R) );
            
        end
        
        
        
       function obj=times(L,R)
            
          obj = RobustSparse( times@double(L,R) );
            
        end

        function obj=mtimes(L,R)
            
          obj = RobustSparse( mtimes@double(L,R) );
            
        end
        
         function obj=rdivide(L,R)
            
          obj = RobustSparse( rdivide@double(L,R) );
            
        end

        function obj=ldivide(L,R)
            
          obj = RobustSparse( ldivide@double(L,R) );
            
        end
        
        function obj=mrdivide(L,R)
            
          obj = RobustSparse( mrdivide@double(L,R) );
            
        end

        function obj=mldivide(L,R)
            
          obj = RobustSparse( mldivide@double(L,R) );
            
        end
   


        
        function obj=inv(L)

          obj = RobustSparse( inv@double(L) );

        end

        function obj=sum(L,varargin)

          obj = RobustSparse( sum@double(L,varargin{:}) );

        end


        function display(obj)
            
            l=inputname(1);
            if isempty(l), l='ans'; end
            disp ' '
            disp([l ' = ']);
            disp ' '
            disp(double(obj));
            
        end
    end
    
end

"Matt " <xys@whatever.com> wrote in message <heeoar$798$1@fred.mathworks.com>...

>
> If you're willing, the following revised tests should give the correct results. Otherwise, it will default with a warning to the normal sparse indexing behavior
>

I sure will, but I cannot take a closer look before Wednesday for many reasons. I'm sure you will make the indexing code right Matt as you seems to be very keen on oop.

Thanks,

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <heepeh$hq1$1@fred.mathworks.com>...

>I'm sure you will make the indexing code right Matt as you seems to be very keen on oop.
====

Well this seems like a particularly good application of it, although hopefully TMW will fix the bug soon and render all this unnecessary.

On 23 Nov, 18:02, "Matt " <x...@whatever.com> wrote:
> "Bruno Luong" <b.lu...@fogale.findmycountry> wrote in message <heefss$lc...@fred.mathworks.com>...
> > %%%%%%%%%%%%%%%%%%
>
> > This is one of the reason why I wrote this package
Thank you for the help!

I also didn't realize that problem...



> >http://www.mathworks.com/matlabcentral/fileexchange/23488-sparse-sub-...
>
> > >> [i j]=find(A);
> > >> AB = setsparse(A,i,j,getsparse(B,i,j),@(a,b) a./b)
>
> > AB =
>
> >              (100000,100000)                  1
>
> > % Bruno
>
> I didn't realize the extent of the problem. It wonder if one could provide a subclass of @double to interface to this package. That way, we could use the old indexing syntax instead of things like this
>
> AB = setsparse(A,i,j,getsparse(B,i,j),@(a,b) a./b)

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <heepeh$hq1$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <heeoar$798$1@fred.mathworks.com>...
>
> >
> > If you're willing, the following revised tests should give the correct results. Otherwise, it will default with a warning to the normal sparse indexing behavior
> >
>
> I sure will, but I cannot take a closer look before Wednesday for many reasons. I'm sure you will make the indexing code right Matt as you seems to be very keen on oop.
======





I carried on with it, and it seems to be doing the job.



>> A=sparse(100000,100000,1), B=RobustSparseDouble(A),

A =

             (100000,100000) 1


B=

             (100000,100000) 1
 



>> whos A B
  Name Size Bytes Class Attributes

  A 100000x100000 400136 double sparse
  B 100000x100000 400084 RobustSparseDouble




>> A(logical(A))=2, %Wrong

A =

              (65408,14101) 2
             (100000,100000) 1




>> B(logical(B))=2, %RIGHT!!!

B=

             (100000,100000) 2
 


I also took care of that pesky problem that native MATLAB sparse matrices have of not being able to interact with singles and other data types:


  %ANNOYING....
>> y=A*single(rand(100000,1)); whos y
  ??? Error using ==> mtimes
Sparse single arrays are not supported.


%BETTER!!!!!!!
>> y=B*single(rand(100000,1)); whos y
  Name Size Bytes Class Attributes

  y 100000x1 800000 double



I guess I'll do seem clean-up then, and put it on the FEX.

Great!!! That's will make the indexing much more intuitive. Thanks.

Bruno

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hgh3t7$c3m$1@fred.mathworks.com>...
       

>
> >> A(logical(A))=2, %Wrong
>
> A =
>
> (65408,14101) 2
> (100000,100000) 1
>

This bug still holds in 2010A prerelease. ;-(

Bruno

RobustSparse works well so far:

>> B(end) % RobustSparse
 
ans =
 
   (1,1) 1

>> A(end) % Matlab sparse
??? Maximum variable size allowed by the program is exceeded.
 
>> >> B(end+1)=1 % RobustSparse
 
B =
 
             (100000,100000) 1
                  (1,100001) 1

>> A(end+1)=1 % Matlab sparse
??? In an assignment A(I) = B, a matrix A cannot be resized.

Bruno

Matt, calling setspase/getsparse on RosbustSparse does not work.

>> setsparse(B,1e5,1e5,10)
??? Error using ==> setspvalmex
First input argument must be a sparse matrix.

Error in ==> spsubsasgn at 120
        S = setspvalmex(S, r(:), c(:), v(:));

Error in ==> setsparse at 45
[vout{:}] = spsubsasgn(varargin{:});
 
>> setsparse(A,1e5,1e5,10)

ans =

             (100000,100000) 10

This is because Matlab API mxIsSparse() return 0 for RobustSparse object. We can make it works, but it requires a change either on your side or mine. Do you accommodate from your side or I do from my side?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hgiapj$p1c$1@fred.mathworks.com>...

> This is because Matlab API mxIsSparse() return 0 for RobustSparse object. We can make it works, but it requires a change either on your side or mine. Do you accommodate from your side or I do from my side?
>
> Bruno

No problem, I can add getsparse/setsparse as RobustClass methods.

An unfortunate difficulty with this "wrapper class" solution is that many superclass functions operating on the RobustSparse subclass objects will unfortunately return ordinary superclass sparse objects, unless I overload those functions, and it's very easy to overlook this when coding. For example:

>> B=RobustSparseDouble( speye(3) ), whos B

B=

   (1,1) 1
   (2,2) 1
   (3,3) 1
 
  Name Size Bytes Class Attributes

  B 3x3 108 RobustSparseDouble

>> BB=triu(B,1), whos BB

BB =

   All zero sparse: 3-by-3

  Name Size Bytes Class Attributes

  BB 3x3 28 double sparse



I've made an effort to overload some of the common functions (sum, inv, plus, minus, etc...) , but it'll be a pain to do them all, and I certainly couldn't keep up with TMW's release of new functions.

This means that one has to be very careful to preconvert a sparse matrix to robust, before doing logical/linear indexing ....[sigh]

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hgh3t7$c3m$1@fred.mathworks.com>...
>
>
> I guess I'll do seem clean-up then, and put it on the FEX.

It's up there now.

http://www.mathworks.com/matlabcentral/fileexchange/26181

Guess we'll see if it's useful to anyone.