function [xstar,ystar,K,y,Sdiag,A,b,t,xi]=phillipsmod(m,n,errorlevel)

% function [xstar,ystar,K,y,Sdiag,A,b,t,xi]=phillipsmod(m,n,errorlevel)
% This function computes the data for the modified Phillips problem
% defined in:
%
%    Residual Periodograms for choosing regularization parameters
%    for ill-posed problems
%    Bert W. Rust and Dianne P. O'Leary
%
% The continuous problem is
%  
%    y(t) = integral_{-3}^3 K(t,xi) x(xi) d xi,  for -3 <= t <= 3,
%
% where
%
%    K(t,xi) = (1 + cos((xi-t) pi/3) )/6 for |xi-t| <= 3 and |t| <= 6,
%            =    0                      otherwise,
% and
%
%    x(xi) = beta(xi) + c1(xi) + c2(xi) + c3(xi),
%
% with
%
%    beta(xi) = alpha0 * (1 + cos(xi pi/3) ) for |xi| <= 3,
%             =   0                          otherwise,
%
% and for k=1,2,3,
%
%    ck(xi) = alpha(k) (1 + cos(2 pi (xi - phi(k)) )) for |xi - phi(k)| <= .5,
%           =   0                                     otherwise.
% 
% The constants alpha0, alpha(k), and phi(k) are defined below.
%
% Input:
%    m = number of observations y
%        (number of    rows in the kernel matrix).
%    n = number of unknown values of x to be determined
%        (number of columns in the kernel matrix).
%    errorlevel defines the size of the random perturbations
%        added to the right-hand side of the problem.
%        (In the paper, errorlevel = 1.e-5.)
% Output:
%    xstar = true solution (n x 1).
%    ystar = "true" right-hand side, computed as Kstar * xstar.
%    K     = discretization of the kernel using trapezoidal
%                rule integration (m x n).
%    y     = ystar + random perturbations with mean 0 and
%            standard deviation S.
%    Sdiag = vector of diagonal elements of S (m x 1),
%            with entries errorlevel * sqrt(ystar).
%    A     = S^{-1} K.
%    b     = S^{-1} y.
%    t     = points at which the data y was measured (m x 1).
%    xi    = points at which the solution y is desired (n x 1).
%
% The resulting problem A x = b is now scaled so that the
% errors in the right-hand side are normally distributed
% with mean zero and standard deviation one.
%
% DPO 07-2007

% Set some constants.

alpha0 = 0.1;
alpha = [.5 .5 1]';
phi = [-1.5 .5 1.5]';
alimit = 5.9625;

% Set the points t and xi.

t   = linspace(-alimit,alimit,m);
xi =  linspace(-3,     3,     n)';

% Set the quadrature weights for trapezoidal rule.

hquad = 6/(n-1);
w = hquad*ones(n,1);
w(1) = w(1)/2; 
w(end) = w(end)/2;

% Set the true solution xstar.

beta = alpha0*(1 + cos(pi*xi/3));
beta = beta .* (abs(xi) <= 3);

for k=1:3,
    c(:,k) = (1 + cos(2*pi*(xi-phi(k))));
    c(:,k) = c(:,k) .* (abs(xi-phi(k)) <= .5);
end

xstar = beta + c * alpha;

% Compute the kernel matrix K.
% (It would be more efficient to do this columnwise.)

K = zeros(m,n);

for i=1:m,
    Kval = (1 + cos((xi-t(i))*pi/3))/6;
    Kval = Kval .* (abs(xi-t(i)) <= 3);
    K(i,:) = Kval .* w;
end

% Compute the right-hand side ystar.

ystar = K * xstar;

% Compute S.
% For values of ystar that are zero, the corresponding
% S(i) is set to machine-epsilon.

Sdiag = errorlevel * sqrt(ystar);
i = (ystar == 0);
if (length(i)>0)
   Sdiag(i) = eps;
end

% Add errors to ystar to obtain y.

y = ystar + Sdiag .* randn(m,1);

% Compute the matrix and right-hand side for the scaled problem.

Sinv = spdiags(1./Sdiag,0,m,m);
A = Sinv * K;
b = Sinv * y;
