function parabolicmeasurecontrol
% This function solves the parabolic measure space control problem
%  min 1/2 |y-yd|_L^2(Q) + alpha |u|_L^2(0,T;M) s.t. y_t - \nu\Delta y = u
% using a semi-smooth Newton method described in the paper
%    "Approximation of elliptic control problems in measure spaces with
%     sparse solutions"
% by Eduardo Casas, Christian Clason and Karl Kunisch, see
% http://math.uni-graz.at/mobis/publications/SFB-Report-2012-007.pdf.
%
% March 28, 2012            Christian Clason <christian.clason@uni-graz.at>
%                                            http://www.uni-graz.at/~clason

%% setup
Nh    = 64;      % number of spatial grid points
Nt    = 256;     % number of time steps
alpha = 1e-1;    % constant blocks of Newton matrix
nu    = 1e-1;    % diffusion coefficient
maxit = 20;      % max iterations in semismooth Newton method

%  setup grid
x = linspace(-1,1,Nh);
t = linspace(0,2,Nt);
[xx,tt] = ndgrid(x,t);
tplot = @(n,f,s) tplot_(n,f,s,Nh,Nt,xx,tt);

% target
yd = tt.*(1-abs(xx));
tplot(1,yd,'target');

% setup finite element discretization
h   = x(2)-x(1);                              % mesh size
tau = t(2)-t(1);                              % time step size
ex = ones(Nh,1);
et = ones(Nt,1);
D2 = spdiags([-ex 2*ex -ex]/h,-1:1,Nh,Nh);    % A_h
Mx = spdiags([ex 4*ex ex]*(h/6),-1:1,Nh,Nh);  % M_h
Dt = spdiags([-et et]/tau,-1:0,Nt,Nt);
It = speye(Nt,Nt);

% space-time discretization
L  = kron(Dt,Mx) + nu*kron(It,D2);  % parabolic operator
Ms = kron(It,Mx);         % mass matrix M_\sigma
II = kron(It,ex);         % summation over spatial coeffs at each time step

%% initialization
p = zeros(Nh*Nt,1);
y = zeros(Nh*Nt,1);
c = alpha*ones(Nt,1);
l = 1;
cq = II*c;                % c(t) as c(t,x)
Ap = (p > cq);    Ap_old = Ap;
Am = (p < -cq);   Am_old = Am;

% constant blocks of Newton matrix
Lyy = -Ms;
Lpy = L;
Lpl = spalloc(Nh*Nt,1,0); Lyl = Lpl;
Lyc = spalloc(Nh*Nt,Nt,0);

% continuation strategy
for gamma = 10.^(0:12)
    disp(['Solving for gamma = ' num2str(gamma,'%1.0e')]);
    
    % semismooth Newton iteration
    for iter = 1:maxit        
        % gradient
        Ly = L'*p - Ms*(y-yd(:));
        Lp = L*y + gamma*(Ap.*(p-cq)+ Am.*(p+cq));
        Lc = 2*l*c + gamma*II'*(Am.*(p+cq) - Ap.*(p-cq));
        Ll = tau*sum(c.^2)-alpha^2;
        
        F = -[Ly; Lp; Lc; Ll];
        
        % Newton matrix
        Lpp = gamma*spdiags(Ap+Am,0,Nh*Nt,Nh*Nt);
        Lpc = gamma*spdiags(Am-Ap,0,Nh*Nt,Nh*Nt)*II;
        Lcc = spdiags(2*l + gamma*II'*(Ap+Am),0,Nt,Nt);
        Lcl = 2*c;
        
        H = [ Lyy Lpy' Lyc      Lyl; ...
            Lpy   Lpp  Lpc      Lpl; ...
            Lyc'  Lpc' Lcc      Lcl; ...
            Lyl'  Lpl' tau*Lcl' 0];
        
        % update
        dx = H\F;
        
        y = y + dx(1:Nh*Nt);
        p = p + dx(1+Nh*Nt:2*Nh*Nt);
        c = c + dx(1+2*Nh*Nt:end-1); cq = II*c;
        l = l + dx(end);
        
        % compute active sets
        Ap = (p > cq);
        Am = (p < -cq);
        
        update = nnz((Ap-Ap_old)|(Am-Am_old));
        ngrad = norm(F(:));
        cviol = norm(max(0,p-cq)+min(0,p+cq));
        fprintf(['Iteration %d: %d points changed in active set, '...
            '|c|^2 = %1.2e, |g| = %e, feasibility: %e\n'], ...
            iter, update, tau*sum(c.^2), ngrad, cviol);
        
        if update == 0
            break
        end
        
        Ap_old = Ap;
        Am_old = Am;
    end
    
    % terminate continuation if Newton iteration converged in one step
    if iter==1
        break;
    end
end

% compute starting point for control
u  = -gamma*(max(0,p-cq)+min(0,p+cq));
Ap = (-u + p - cq > 0);   Ap_old = Ap;
Am = (-u + p + cq < 0);   Am_old = Am;

%% solve problem
disp('Compute optimal control');

% constant blocks of Newton matrix (if changed from above)
Luu = -speye(Nh*Nt,Nh*Nt);
Luy = L;
Lup = spalloc(Nh*Nt,Nh*Nt,0);  Lyu = Lup;  Lpy = Lup;
Luc = spalloc(Nh*Nt,Nt,0);     Lyc = Luc;
Lul = spalloc(Nh*Nt,1,0);      Lyl = Lul;  Lpl = Lul;
Lyy = -Ms;

% semismooth Newton iteration
for iter = 1:maxit    
    % gradient
    Lu = L*y - u;
    Ly = L'*p - Ms*(y-yd(:));
    Lp = u + (Ap.*(-u+p-cq) + Am.*(-u+p+cq));
    Lc = 2*l*c + II'*(Am.*(-u+p+cq) - Ap.*(-u+p-cq));
    Ll = tau*sum(c.^2)-alpha^2;
    
    F = -[Lu; Ly; Lp; Lc; Ll];
    
    % Newton matrix
    Lpp = spdiags(Ap+Am,0,Nh*Nt,Nh*Nt);
    Lpu = -Luu-Lpp;
    Lpc = spdiags(Am-Ap,0,Nh*Nt,Nh*Nt)*II;
    Lcu = -Lpc;
    Lcc = spdiags(2*l + II'*(Ap+Am),0,Nt,Nt);
    Lcl = 2*c;
    
    H = [ Luu  Luy  Lup  Luc      Lul;...
        Lyu    Lyy  Luy' Lyc      Lyl; ...
        Lpu    Lpy  Lpp  Lpc      Lpl; ...
        Lcu'   Lyc' Lpc' Lcc      Lcl; ...
        Lul'   Lyl' Lpl' tau*Lcl' 0];
    
    % update
    dx = H\F;
    
    u = u + dx(1:Nh*Nt);
    y = y + dx(1+Nh*Nt:2*Nh*Nt);
    p = p + dx(1+2*Nh*Nt:3*Nh*Nt);
    c = c + dx(1+3*Nh*Nt:end-1);  cq = II*c;
    l = l + dx(end);
    
    % compute active sets
    Ap = (-u + p - cq > 0);
    Am = (-u + p + cq < 0);
    
    update = nnz((Ap-Ap_old)|(Am-Am_old));
    ngrad  = norm(F(:));
    fprintf(['Iteration %d: %d points changed in active set, '...
        '|c|^2 = %1.2e, |g| = %e\n'], ...
        iter, update, tau*sum(c.^2), ngrad);
    if update == 0
        break
    end
    
    Ap_old = Ap;
    Am_old = Am;
end

%% plot results
tplot(2,u,'control');
figure(3);plot(t,c);xlabel('t');ylabel('c_\sigma(t)');title('norm bound');
tplot(4,y,'state');


function tplot_(n,f,s,Nh,Nt,xx,tt)
figure(n); surf(tt,xx,reshape(f,Nh,Nt)); axis tight
shading interp;colormap summer;lighting phong;camlight headlight;alpha(0.8)
title(s); xlabel('t'); ylabel('x');
drawnow;
