function [e hatA ek]=A_welasso(RT60,Aind,X,Sinf,pena)
%performs several FISTA algorithms to solve constrained problem (3)
%							\min_\A  p(\A)
%							s.t  \|\X-\A\star \s\|_2 = 0
%
%
%Input
%		X		M*T+K-1		Observations matrix
%		Sinf	N*T			Sources matrix
%		Aind	M*N*K		The original filter, only used for the peak value \tD of the expected time enveloppe \rho
%		RT60				Reverberation time (ms) of the filter
%		pena				integer between 1 and 5 for the choice of the penalisation (pena=1 for P_{1,\rho}, pena=2 for P_{2,\rho}, pena=3 for P_1, pena=4 for P_2, pena=5 for P_{1,2})
%                  
%
%Output
%		hatA	M*N*K		the estimation of the filter
%		hatAs	M*N*K*14	estimated filters after each of the 14 iterations of FISTA


options.maxit = 500;


[M,useless]=size(X);
[M,N,K]=size(Aind);

options.Aorig=Aind;

% choix de la liste de lambda
lambda_vect = [1e-15,1e-14,1e-13,1e-12,1e-11,1e-10,1e-9,1e-8,1e-7,1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1];
nbxp = numel(lambda_vect);
lambda_vect=lambda_vect(end:-1:1);

% diverses initialisations
hatAs=zeros(M,N,K,nbxp);

X(:,end-K+1:end)=0;
%tD est l'abscisse du plus grand pic de A pour initialiser l'enveloppe temporel \rho (initialisation oracle ici)
[val,tD]=max(abs(Aind),[],3);
tD=floor(mean(tD(:)));

fs= 44100 ; % sampling frequency

noscale=ones(M,N,K);
scaleamp=ones(M,N,K);


RT60c=RT60;
options.itplot=[];
for n=1:N
    for m=1:M
        [val,tD]=max(abs(Aind(m,n,:)),[],3);
        %v_dp=max(abs(Aind(m,n,tD+3:end)));
        tD=tD+3;
        tD=1;
        namp=-(0:K-tD-1)/fs*1000/RT60c*60; % amplitude in dB
        namp=[10.^(namp/20)]; % amplitude in linear scale
        
        scaleamp(m,n,tD+1:K)=namp;
    end
end

switch pena
    case 1
        options.tR=K-1;%L1
        options.rescale=scaleamp;%rescale on
        
    case 2
        options.tR=1;%L2
        options.rescale=scaleamp;%rescale on
        
    case 3
        options.tR=K-1;%L1
        options.rescale=noscale;%rescale off
        
    case 4
        options.tR=1;%L2
        options.rescale=noscale;%rescale off
        
    case 5
        
        options.tR=floor(K/3); %60ms after the peak tD
        options.rescale=scaleamp;%rescale on
    case 6
        options.rescale=noscale;%rescale off, won't be used
        lambda_vect=0*lambda_vect;
end


k=0;
for ll =  lambda_vect
    k=k+1;
    options.lambda=ll;
    
    [hatA L]=inverse_A_time(X,Sinf,options);
    options.init=hatA;
    options.lipschitz=L;
    ek(k)=SRRA(Ahp(hatA),Ahp(Aind));
end
e=ek(end);
end


% Written by Matthieu Kowalski, 2010
% Updated by Alexis Benichoux, 2013