function gf=gradloglikGaP(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model 
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt


Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

% % % Hkt_linear=x(1:end-peval.ncomp*2); % intensities
% % % cx=x(end-peval.ncomp*2+1:end-peval.ncomp); %x-coordinates of the centers
% % % cy=x(end-peval.ncomp+1:end); % y-coordinates of the centers
% % %
Hkt_linear=x; % intensities
cx=varargin{6};
cy=varargin{7};
% % %

sigpsf_vec=repmat(sigpsf,peval.ncomp,1); %all psfs same sigma
cxy_vec=[cx'+.5, cy'+.5];
a_vec=1./(sigpsf_vec*2*pi); % all normalised to 1

Hkt=reshape(Hkt_linear, peval.ncomp, peval.nt);
% generate PSFs from given parameters:
Wxkpix=gauss2dmultislice([peval.nx, peval.ny, peval.ncomp], cxy_vec, sigpsf_vec, a_vec);
Wxkpix=normalizePSF(Wxkpix); %normalize PSFs to 1
Wxk=reshape(Wxkpix,peval.nx*peval.ny, peval.ncomp);
% add the background component (just for muttiplication - in the gradient it should not matter): 
[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
P=Wxkbg*Hktbg; %current approximation

xxvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cx, 'xx'); %linera grasdient shifted by cx
yyvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cy, 'yy');

% dW/dcx:
Wxtcx=1/sigpsf^2*xxvc.*Wxk; 
% dW/dcy:
Wxtcy=1/sigpsf^2*yyvc.*Wxk;

% d(log(L))/dHkt:
% % % gfHkt=(alpha-1)*1./Hkt - 1/beta + Wxk'*(Vxt./P-eye(peval.nx*peval.ny,peval.nt));
% % % gfHkt= Wxk'*(Vxt./P-eye(peval.nx*peval.ny,peval.nt));
gfHkt= Wxk'*(Vxt./P)-1; %without background
% d(log(L))/dcx:
gfcx=diag(Wxtcx'*(Vxt./P)*Hkt'); 
% d(log(L))/dcy:
gfcy=diag(Wxtcy'*(Vxt./P)*Hkt'); 

% % % gf = [reshape(gfHkt',1,peval.nt*peval.ncomp), gfcx', gfcy'];
% gf = [reshape(gfHkt(1:peval.ncomp, :),1,peval.nt*peval.ncomp)];
gf = [reshape(gfHkt,1,peval.nt*peval.ncomp)];

end
function Wnorm=normalizePSF(W)
sw=size(W);
Wr=reshape(W, sw(1)*sw(2),sw(3));
q=squeeze(sum(Wr,1));
Wrnorm=Wr./repmat(q,sw(1)*sw(2),1);
Wnorm=reshape(Wrnorm,sw(1), sw(2), sw(3));
end
function xxvc = lineargrad(sizevec, cx, dir)
switch dir
    case 'xx'
        xxp=double(xx(sizevec, 'true')); %linear function - pixels
    case 'yy'
        xxp=double(yy(sizevec, 'true')); %linear function - pixels
    otherwise 
        error('Wrong dir')        
end
    
xxv=reshape(xxp,sizevec(1)*sizevec(2),sizevec(3)); %linear function - vector
xxvc=xxv-repmat(cx,sizevec(1)*sizevec(2),1);
% xxp=double(xx([peval.nx, peval.ny, peval.ncomp], 'true')); %linear function - pixels
%yyp=double(yy([peval.nx, peval.ny, peval.ncomp], 'true'));
% xxv=reshape(xxp,peval.nx*peval.ny,peval.ncomp); %linear function - vector
%yyv=reshape(yyp,peval.nx*peval.ny,peval.ncomp);
% xxvc=xxv-repmat(cx,peval.ncomp,1);
%yyvc=yyv-repmat(cy,peval.ncomp,1);
end