function [xopt,pen_final] = commGSP_quad_regularization(A, sig, y_sig, pens, type)
% Optimizes over x the objective ||y_sig-x||_2^2 + pen .* x*P*x^T
%   (may be used for signal denoising under the assumption of
%   modularity, smoothness, etc. depending on P)
% INPUT:
%   A - (symmetric) adjacency matrix
%   sig - ground truth signal (if uknown, put the same as y_sig)
%   y_sig - input noisy signal to the optimization/denoising
%   pens - vector of penalty values to explore; the procedure
%          automatically searches one level deeper meaning that it checks
%          length(pens) more values interpolated around the value of
%          pens found as optimal during the first search
%   type - string defining matrix P in the 
%          regularization term x*P*x^T:
%          'laplacian' - puts L (smooth signal reconstruction)
%          'modular' - puts spectrally shifted modularity
%                      (modular signal reconstruction)
%          'anti-modular' - puts spectrally shifted modularity
%                           (anti-modular signal reconstruction)

% OUTPUT
%   xopt - reconstructed/denoised signal
%   pen_final - chosen (optimal) value of the penalty parameter

N=size(A,1);
if strcmp(type, 'laplacian')
    P = commGSP_build_matrix(A, 'laplacian');
else
    Q = commGSP_build_matrix(A, 'modularity');
    [~, LamQ] = commGSP_eig(Q, 'descend');
    if strcmp(type, 'modular')
        P = -Q + eye(N)*max(diag(LamQ));
    end
    if strcmp(type, 'anti-modular')
        P = Q - eye(N)*min(diag(LamQ));
    end
end

options = optimoptions('quadprog');
options.Display = 'off';

N=size(P,2);
tmp_L=sparse(N,length(pens));
for pens_iter=1:length(pens)
    tmp_L(:,pens_iter) = quadprog(2.*pens(pens_iter).*P+eye(N), -2*y_sig,[],[],[],[],[],[],[],options);
    tmp_fval_L(pens_iter)=sqrt(mean((sig-tmp_L(:,pens_iter)).^2));
end
[~,ii]=min(tmp_fval_L);

if ii==1; pen_low=1; else; pen_low=ii-1; end
if ii==length(pens); pen_high=length(pens); else; pen_high=ii+1; end

pens_new=linspace( pens(pen_low), pens(pen_high), length(pens));
for pens_iter=1:length(pens)
    tmp_L(:,pens_iter) = quadprog(2.*pens_new(pens_iter).*P+eye(N), -2*y_sig,[],[],[],[],[],[],[],options);
    tmp_fval_L(pens_iter)=sqrt(mean((sig-tmp_L(:,pens_iter)).^2));
end

[~,ii]=min(tmp_fval_L);

xopt = tmp_L(:,ii);
pen_final = pens_new(ii);

end