function [ xc, Nc, resn ] = ALqMin( fw, U, q, tol, No, itmax )
% ALqMin adaptive lq like minimization
%
% It gives an approximated solution to an underdetermined least square problem by
% minimzation of the lq^q-like quatity. The algorithm evolves by adaptive selection of a
% subset of normal equations as linear constraints for the minimization process [1].  The
% constrained minimization is implemented by the function LQ_FOCUSS, which applies the
% FOCUSS algorithm [2].
%
% Usage:    [ xc ] = ALqMin( fw, U, q );
%           [ xc, ki, resn ] = ALqMin( fw, U, q, tol, No, itmax );
%
% Inputs:
%   fw      data to be modeled  by fw~U*xc 
%   U       matrix in the model above
%   q       specifies the value q for the lq norm
%   tol     the routine will stop if norm(fw-U*xc)<tol  (default tol=1e-8)
%   No      maximum number of constraints (stopping condition) (default No=size(U,2)) 
%   itmax   maximal number of iterations for the focuss algorithsm (default itmax=30) 
%
% Outputs:
%   xc      solution of the lq minimization
%   Nc      total number of normal equations that have been considered
%   resn    value of norm(fw-U*xc)*sqrt(delta)
%
% References:
%  [1] L. Rebollo-Neira and A. Plastino, Nonlinear non-extensive approach for
%  identification of structured information,  Physica A, 2009
%  [2] B.D. Rao and K. Kreutz-Delgado,  An Affine Scaling Methodology for Best Basis 
%  Selection, IEEE Trans. Sig. Proc. 47 (1999) 187--200

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/

  if nargin<6,   itmax=30; end
  if(nargin<5) | (isempty(No)==1); No=size(U,2); end
  if(nargin<4) | (isempty(tol)==1); tol=1e-8; end

name='ALQMIN';
[L,N]=size(U);
G=U'*U;
mu=U'*fw';
clear fp x0
CA=sum(fw*U)/sum(sum(G));
for n=1:N; x0(n)=CA;end %initializing the iterative precess
fp=U*x0';
mup=U'*fp;
errmu=abs(mu-mup);
[max_e,ind]=max(errmu);
H=min(No,N);
phi=[];
indp=[];
muk=[];
itmax=50;
tolf=1e-7;
for k=1:H;
indp(k)=ind;
%
%adds one more constraint
%
muk(k)= mu(ind);
phi(k,:)=G(ind,:);
x0=pinv(phi)*muk';
%
%solves by focuss
%
[xc]=LqFOCUSS(muk',phi,x0,q,tolf,itmax);
fp=U*xc;
mup=U'*fp;
nop(k)=sum(abs(xc).^q);
%
resn=norm(fp'-fw); %norm of the errors between model and data
%
if (resn <tol) break; end;
mup=U'*fp;
errmu=abs(mu-mup);
[max_e,ind]=max(errmu);%choses new constraint
end
Nc=k;%total number of constraints 

%Laura Rebollo-Neira (2009)
