%A temporary replacement function until bsxfun supports sparse matrices natively (see bug #41441)
%sig is logical 3x3 matrix that indicates how the function op behaves on negative, zero, and positive values
%(where the index 1 is for negative, 2 is for zero, and 3 is for positive). sig(i,j) should be true if and only if
%the value of op can be nonzero. eg if sig(1,2) is true, that means there exists some x, y such that x < 0 and y == 0
%and op(x,y) ~= 0
%It is assumed that op behaves properly on matrices of the same size and also on a matrix together with a singleton
function [res] = spbsxfun(op,A,B,sig)
if size(A) == size(B) || size(A) == [1,1] || size(B) == [1,1]
res = op(A,B);
return;
endif
msize = max(size(A),size(B));
% TODO Change asserts to dimension mismatch exception
assert(size(A) == msize | size(A) == [1,1]);
assert(size(B) == msize | size(B) == [1,1]);
res = sparse(msize(1),msize(2));
for i = 1:3
for j = 1:3
if sig(i,j)
rel = relevant(A,compfuncs{i},B,compfuncs{j});
arel = spbsxmult(A,rel);
brel = spbsxmult(B,rel);
res += op(arel,brel);
%TODO These last three lines only work if op(0,0) == 0
%I would like to change them to:
%res(rel) += op(A(rel),B(rel));
%except that logical indexing into sparse matrices is
%broken (see bug #40341)
endif
endfor
endfor
%TODO Convert res to the correct return class for op
endfunction
compfuncs = {@lt,@eq,@gt};
%Finds the set of all entries for which cfA(A,0) & cfB(B,0) holds.
%This is optimized to look first in whichever of A or B is a vector
%and then only examine the relevant values in the matrix argument
function [res] = relevant(A,cfA,B,cfB)
ismatA = min(size(A)) > 1;
ismatB = min(size(B)) > 1;
if ismatA
assert(~ismatB);
mat = A;
cfmat = cfA;
other = cfB(B,0);
elseif ismatB
assert(~ismatA);
mat = B;
cfmat = cfB;
other = cfA(A,0);
else
relA = cfA(A,0);
relB = cfB(B,0);
if(size(A,1) == 1)
assert(size(B,2) == 1);
row = relA;
col = relB;
else
assert(size(A,2) == 1);
assert(size(B,1) == 1);
row = relB;
col = relA;
endif
res = col * row;
return;
endif
res = sparse(size(mat,1),size(mat,2));
if size(other,1) == 1
relmat = mat(:,other);
relmat = cfmat(relmat,0);
res(:,other) = relmat;
else
assert(size(other,2) == 1);
relmat = mat(other,:);
relmat = cfmat(relmat,0);
res(other,:) = relmat;
endif
endfunction
function [res] = spbsxfun(A,logicmat)
if(size(A,1) == 1)
%TODO
%Using spdiag rather than diag in case A is logical
%since diagonal logical matrices are not supported.
%See bug #40105
res = logicmat * spdiag(A);
elseif(size(A,2) == 1)
%TODO
%Using spdiag rather than diag in case A is logical
%since diagonal logical matrices are not supported.
%See bug #40105
res = spdiag(A) * logicmat;
else
res = A .* logicmat;
endif
endfunction
function [r] = spdiag(v)
assert(any(size(v) == 1));
v = reshape(v,[],1);
r = spdiags(v,0,length(v),length(v));
endfunction