function F = predictor_matrix(X, thetas, f, regtype, lambda)
% F = predictor_matrix(X, thetas, f, regtype, lambda)
%  Generates the ? x K matrix of regressors F, (either n x K or (n+1) x K)
%  
% Inputs:
%  - X, p x n design matrix X
%  - thetas, D x K matrix of parameters
%  - f, function handle of the form f(x,theta): specifies kernel family
%  - regtype, 'none', 'L1', or 'unweighted': regularization type
%  - lambda, scalar: penalization constant
%
% Outputs:
%  - F, ? x K matrix.  A n x K matrix if regtype == 'none', 
%     a (n+1) x K matrix if regtype == 'unweighted' or 'L1'.
%     Signal matrix with F(i,j) = f(X(:,i), theta(:,j))
%
% Example usage:
%  theta_space = .5 .* (fullfact([5 5])-3)'; % grid on [-1,1]^2
%  delta = .5;
%  thetas = sample_parameters(theta_space,delta,10); % will sample points from [-1.25,1.25]^2
%  f = @(x,theta) exp(-norm(x-theta)^2); % gaussian kernel
%  X = randn(2,10);
%  regtype = 'none'; lambda = 0;
%  F = predictor_matrix(X,thetas,f,regtype,lambda);

K = size(thetas,2);
n = size(X,2);
rawF = zeros(n,K);
fhandle = @(i,j) f(X(:,i),thetas(:,j));
cls = ones(n,1) * (1:K);
rws = (1:n)' * ones(1,K);
rawF = arrayfun(fhandle, rws, cls);

if strcmp(regtype,'none');
  F = rawF;
end

if strcmp(regtype,'L1') || strcmp(regtype,'unweighted');
  F = [rawF; sqrt(lambda).*ones(1,K)];
end

end
