function losses=testWs(X,y,wrange,lossfunc)
% losses=testWs(X,y,wrange,lossfunc)
% Arguments:
%  X - N x d+1 design matrix
%  y - observed function values
%  wrange - a cell array, with d+1 cells. wrange{j} is the range over
%           which the j-th coefficient has to run
%  lossfunc - name of the function used to compute empirical loss.
%             lossfunc is assumed to take two arguments: yhat (estimated)
%             and y (observed).
% 
% Output:
%  losses(i,j,...) is the value of the loss for the corresponding set of
%  parameter (regression coefficient) values

if (size(X,2) ~= length(wrange))
  error('Sizes of X and wrange should match');
end

% this is only working for up to 3D
if (length(wrange)==2)
  [w0,w1]=meshgrid(wrange{1},wrange{2});
  ws=[reshape(w0,1,[]);reshape(w1,1,[])];
else
  [w0,w1,w2]=meshgrid(wrange{1},wrange{2},wrange{3});
  % order is switched - see help meshgrid for explanation.
  ws=[reshape(w1,1,[]);reshape(w0,1,[]);reshape(w2,1,[])];
end

yhat=X*ws;

loss=zeros(1,size(yhat,2));
% NOTE: this loop could be avoided by making sure that the loss
% calculating functions know how to deal with matrices, rather than
% vectors, as input. (you could then pass the entire yhat and y)
for i=1:length(loss)
  loss(i)=feval(lossfunc,yhat(:,i),y);
end
if (length(wrange)==2)
  losses=reshape(loss,length(wrange{2}),length(wrange{1}))';
else
  losses=reshape(loss,length(wrange{1}),length(wrange{2}),[]);
end