%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
