function [y,y0] = generate_signal(X, thetas, w, f, noise_type, noise)
% [y,y0] = generate_signal(X, thetas, w, f, noise_type, noise)
%  Generates a noisy signal and noiseless signal from a mixture model, with family f
%
% Inputs:
%  - X, p x n matrix of design points
%  - thetas, 1 x K list of parameters
%  - w, a 1 x K vector of weights
%  - f, function handle of the form f(x,theta): specifies kernel family
%  - noise_type: 'gaussian', 'rician', 'poisson', or 'none', see below:
%  - noise, scalar: scale parameter for generating noise
%    * if noise_type=='gaussian': the standard deviation
%    * if noise_type== 'rician': sqrt(2) * scale parameter, approaches std deviation in the limit to zero
%    * if noise_type== 'poisson': y will be generated as Poisson(y0/noise), y0 will be scaled to y0/noise
% 
% Outputs:
%  - y, n x 1 vector: noisy signal
%  - y0, n x 1 vector: noiseless or 'true' signal
%
% 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
%  w = rand(10,1);
%  f = @(x,theta) exp(-norm(x-theta)^2); % gaussian kernel
%  X = randn(2,10);
%  regtype = 'none'; lambda = 0;
%  noise_type = 'gaussian'; noise=0.1; % specifies gaussian noise N(0, 0.01)
%  [y,y0] = generate_signal(X, thetas, w, f, noise_type, noise);


n = size(X,2);
K = size(w,2);
F = predictor_matrix(X, thetas, f, 'none', 0);
y0 = F * w;
y = y0;
if strcmp(noise_type,'gaussian')
  y = y0 + noise.* randn(n,1);
end
if strcmp(noise_type,'rician')
  y = sqrt((y0+(1/sqrt(2)).*noise.*randn(n,1)).^2 + 1/2.*(noise.*randn(n,1)).^2);
end
if strcmp(noise_type,'poisson')
  y0 = y0./noise;
  y = possrnd(y0,n,1);
end
end
