function  [hatA] = inverse_A_time(X,S,options)
%performs FISTA algorithm to solve (2) 
%							\min_\A\left\{\frac{1}{2}\|\X-\A\star \s\|_2^2+\lambda\p(\A)\right\}	
%Input
%		X M*T+K-1	Observations matrix
%		S N*T		Sources matrix
%		option		Object containing :
%
	%					option.lambda		\lambda constant of problem (2)
	%					option.maxit		number of Fista iterations
	%	(see prox_A)	option.tR			\p=ell_1 norm for t<tR+1 and \p=\ell_2 norm for t>tR
	%	(see prox_A)	option.rescale		\rho(t): expected time enveloppe of filters
	%	(optional)		option.maxitl		number of power_A_iteration iterations for Lipschitz constant estimation
	%	(optional)		option.init			initialization of A for FISTA
	%	(optional)		option.lipschitz	lipschitz constant of \nabla A\mapsto \frac{1}{2}\|\X-\A\star \s\|_2^2
	%	(optional)		option.Aorig		Original filter for perfomance display during iterations
%Output
%
%	    hatA		M*N*K matrix		Estimated filter, Argmin of (2)
%
%

lambda=options.lambda;


[M,TpKm1]=size(X);

[N,T] = size(S);
K=TpKm1-T+1;



if ~isfield(options, 'init')
    options.init = ones(M,N,K);
end
if ~isfield(options, 'maxitl')
	options.maxitl=500;
end
if ~isfield(options, 'lipschitz')
  %% Choice of C: square of the upper bound norm of the operator.
    
      [C Cmin] = power_A_iteration(S,M,K,options.maxitl);
      options.conditionnement=C/Cmin;

  options.lipschitz = C;
end



options


    fista_aux();






function fista_aux()

%% Some initializations

t = 1;
y = options.init;
dccA=y;
iC = 1/options.lipschitz;


%% Main loop: fista iterations 
for it=1:options.maxit
    old = dccA;
    yg=compute_grad_A_time_aux(y);
    dccA = y - iC*yg;
    dccA = prox_A(dccA,options);
    
    oldt = t;
    t = (1+sqrt(1+4*t^2))/2;
    y = dccA + (oldt-1)/t * (dccA-old);
    
 
    if mod(it,floor(options.maxit/2))==0
		if isfield(options, 'Aorig')

			perf=SRRA(options.Aorig,y);  
			disp({it,'SRR_A=',perf});
		else
			disp(it)
		end	
	end

end
hatA=y;
end



function [grad]=compute_grad_A_time_aux(dcc)


    AS =  creation_mel_conv(dcc,S);

    grad=-adjoint_conv(S,X-AS);
end

end
% Written by Matthieu Kowalski, 2010
% Updated by Alexis Benichoux, 2011