function measurecontrol
% This function solves the regularized sparse problem
%  min 1/2 |y-z|_(omObs)^2 + alpha |u| + 1/(2c) |u|^2 s.t. Ay=chi_omCon u
% where omObs and omCon are subdomains, using a semi-smooth Newton method
% and continuation in c, described in the paper
%          "A measure space approach to optimal source placement"
% by Christian Clason and Karl Kunisch, see
% http://math.uni-graz.at/mobis/publications/SFB-Report-2011-015.pdf.
%
% June 22, 2011             Christian Clason <christian.clason@uni-graz.at>
%                                            http://www.uni-graz.at/~clason

%% setup
N     = 128;        % grid size
alpha = 1e-4;       % control penalty
maxit = 10;         % max iterations in semismooth Newton method
cgits = 100;        % max iterations in Krylov method
tol   = 1e-9;       % tolerance for SSN and Krylov method

% coefficients in differential operator
nu = 0.1;           % diffusion
b  = [0,1];         % convection

% setup grid
[xx,yy] = meshgrid(linspace(-1,1,N));
h = 2/N;

% setup subdomains
omObs = ((xx.^2+yy.^2)<1/32);                       % observation domain
omCon = ((xx.^2+yy.^2)>1/16)&((xx.^2+yy.^2)<1/2);   % control domain
niceplot = @(n,f,s) niceplot_(n,f,s,N,xx,yy,omObs,omCon);
omObs = omObs(:); omCon = omCon(:);

% target
z = yy(:).*omObs;

niceplot(1,z.*omObs,'target, control and observation domains');
text(-.05,.0,'\omega_o','FontSize',12);text(.25,-.3,'\omega_c','FontSize',12);

% setup operators
e  = ones(N,1);
D  = spdiags([-e e]/h,-1:0,N,N);
D2 = spdiags([e -2*e e]/h^2,-1:1,N,N);
I  = speye(N,N);
A  = -nu*(kron(I,D2)+kron(D2,I)) - (b(1)*kron(I,D) + b(2)*kron(D,I));
[L,U] = lu(A); Lt = L'; Ut = U';  % precompute LU factorization of A and A'

S  = @(u) omObs.*(U\(L\(omCon.*u)));      % forward operator
Sa = @(p) omCon.*(Lt\(Ut\(omObs.*p)));    % adjoint operator

%% Semismooth Newton iteration
u = zeros(N*N,1);
% continuation strategy
for c=10.^(0:12)
    fprintf('Compute minimizer for c=%1.0e\n',c);
    
    iter   = 0;
    update = 1;
    
    % semismooth Newton iteration
    while (update)&&(iter<maxit)
        
        y = S(u);       % solve state equation
        p = Sa(y-z);    % solve adjoint equation
        
        % compute active sets
        Ap = c*(p > alpha);
        Am = c*(p < -alpha);
        
        F  = u + Ap.*(p - alpha) + Am.*(p + alpha);
        DF = @(x) x + (Ap+Am).*Sa(S(x));
        
        % solve for Newton step
        [du,~] = gmres(DF,-F(:),[],tol,cgits);
        u = u + du;
        
        if iter==0
            update = 1;
        else
            update = nnz((Ap-Ap_old)|(Am-Am_old));
            fprintf(['Iteration %d: %d points changed in active set,\t'...
                'norm of gradient %e\n'], iter, update, norm(F(:)));
        end
        
        iter   = iter+1;
        Ap_old = Ap;
        Am_old = Am;
        
    end
    
end

%% plot results
niceplot(2,u,'control');
niceplot(3,U\(L\(omCon.*u)),'state');

end % main function

function niceplot_(n,f,s,N,x,y,om1,om2)
% plot surfaces with proper shading and indication of subdomains
figure(n); f = reshape(f,N,N);
surf(x,y,f); set(gca,'FontSize',12);
shading interp; camlight headlight; lighting gouraud; alpha(0.67);
hold on;
domainplot(x,y,f,om1,'r');
domainplot(x,y,f,om2,'k');
hold off
title(s); xlabel('x_1'); ylabel('x_2');
drawnow;
end % niceplot_ function

function domainplot(xi,yi,zi,om,linecolor)
% superimposed contour plot of OM on surface plot of ZI
c_xz = contours(xi,zi,om,[1,1]);
c_yz = contours(yi,zi,om,[1,1]);
i = 1;    n = 1;
if isempty(c_xz) || isempty(c_yz)
    return
end
while 1
    lev(n) = c_xz(1,i); %#ok<*AGROW>
    nvals  = c_xz(2,i);
    x{n}   = c_xz(1,i+1:i+nvals);
    y{n}   = c_yz(1,i+1:i+nvals);
    z{n}   = c_xz(2,i+1:i+nvals);
    i = i+nvals+1;
    if i > size(c_xz,2)
        break
    else
        n = n+1;
    end
end
for n=1:length(x)
    plot3(x{n},y{n},z{n},'color',linecolor,'userdata',lev(n),...
        'LineWidth',1.2); hold on
end
end % domainplot function
