function lnq_GD_minimal_alllen(output_name,gen_seq,data_seq,his,max_it,tol,GD_step_aa,theta1)
%likelihood max for 1point model using ME and gradient descent
M=50; %number of exp iterations after which the status is displayed
%nannorm=@(x) sqrt(nansum(x(:).^2));

%check arguments
if ~exist('data_seq','var')
    disp('must supply data pool');
    return;
end;


if ~exist('gen_seq','var')
    disp('must supply gen pool');
    return;
end;


if ~exist('output_name','var') || ~ischar(output_name)
    disp('must supply output name as string');
    return;
end;

if ~exist('his','var')
    his=false;
end;

if ~exist('max_it','var')
    max_it = 100000;
end;

if ~exist('GD_step_aa','var')
    GD_step_aa=.1;
end;

GD_step_L=.1;

if ~exist('tol','var')
    tol=1e-4; %tolerance for change in selection factors suffcient to quit
end;

if iscellstr(gen_seq) %if gen_seq is aa seqs in string format
    gen_seq=uint32(aa2int(char(gen_seq))); %transform seqs to fixed length int arrays
end;
if iscellstr(data_seq)
    data_seq=uint32(aa2int(char(data_seq)));
end;

gen_len=sum(gen_seq~=0,2);
data_len=sum(data_seq~=0,2);
% gen_len=cellfun(@length,gen_seq);
% data_len=cellfun(@length,data_seq);

gen_mask=gen_seq>0; %places of actual values (not padded)
data_mask=data_seq>0;

if max(gen_seq(:))>20
    max_aa=64;
    pd=4;
else
    max_aa=20;
    pd=5;
end;

gen_seq(gen_seq==0)=pd; %padding with unused value instead of 0 so indexing will work
data_seq(data_seq==0)=pd;

gen_samp=size(gen_seq,1);
data_samp=size(data_seq,1);

gen_maxL=size(gen_seq,2);
data_maxL=size(data_seq,2);

maxL=max([gen_maxL data_maxL]);

gen_pos=repmat(uint32(1:gen_maxL),gen_samp,1); %matrix for the posisiton of every aa in every seq
gen_L=repmat(uint32(gen_len),1,gen_maxL); %matrix for the seq length for every aa in every seq
gen_numseq=repmat(uint32(1:gen_samp)',1,gen_maxL); %matrix for the seq # for every aa in every seq
gen_numseq=gen_numseq(gen_mask);

data_pos=repmat(uint32(1:data_maxL),data_samp,1);
data_L=repmat(uint32(data_len),1,data_maxL);

gen_inds={gen_seq(gen_mask),gen_pos(gen_mask),gen_L(gen_mask)}; %3d index (aa, pos, len) for every aa in every seq

data_subs=sub2ind([max_aa maxL maxL],data_seq,data_pos,data_L);
gen_subs=sub2ind([max_aa maxL maxL],gen_seq,gen_pos,gen_L); %linear ind for (aa, pos, len) for every aa in every seq

gen_aa_f=accumarray(gen_inds,1,[max_aa,maxL,maxL])/gen_samp;

%restricting to data seqs that have assig and are represented in gen
data_1ok=all(gen_aa_f(data_subs)>0 | ~data_mask,2); %data seqs having aa/pos/len comb which are in gen

data_ok=data_1ok;

data_seq=data_seq(data_ok,:); %restricting data pool
data_mask=data_mask(data_ok,:);
data_pos=data_pos(data_ok,:);
data_L=data_L(data_ok,:);
data_len=data_len(data_ok);
data_samp=length(data_len);


pc = min(1/data_samp,1/gen_samp); %pseudocount 

data_aa_f = accumarray({data_seq(data_mask),data_pos(data_mask),data_L(data_mask)},1,[max_aa,maxL,maxL])/data_samp + pc;
gen_aa_f = gen_aa_f + pc;

for l=1:maxL;gen_aa_f(:,(l+1):end,l)=0;end;
for l=1:maxL;data_aa_f(:,(l+1):end,l)=0;end;

gen_dlen=hist(gen_len,1:maxL)'/gen_samp;
data_dlen=hist(data_len,1:maxL)'/data_samp;


data_aa_f_N=bsxfun(@rdivide,data_aa_f,shiftdim(data_dlen,-2));
data_aa_f_N(isnan(data_aa_f_N)) = 0;

gen_aa_f_N=bsxfun(@rdivide,gen_aa_f,shiftdim(gen_dlen,-2));
gen_aa_f_N(isnan(gen_aa_f_N)) = 0;

ratio_aa=data_aa_f./gen_aa_f; %enrichment values
ratio_aa_N=data_aa_f_N./gen_aa_f_N; %taking out length dependence from intial value of theta1

ratio_aa(isnan(ratio_aa))=1;
ratio_aa_N(isnan(ratio_aa_N))=1;

Lrange=data_dlen*data_samp>100; %lengths where there is enough data statistics

save([output_name '_parameters']);
%% initilze vars

if ~exist('theta1','var')
    theta1=log(ratio_aa_N);
    theta1(:,:,~Lrange)=0; %set all selection factors in low statistics zone to zero
    thetaL=log(data_dlen./gen_dlen);
end;


if his
    save_it = floor(max_it/M);
    Rtheta1 = zeros([size(theta1) save_it]);
    Rgrad_likeli_aa = zeros([size(theta1) save_it]);
    Rgrad_const_aa = zeros([size(theta1) save_it]);
    Rdelta_aa = zeros([size(theta1) save_it]);
    
    RthetaL = zeros([length(thetaL) save_it]);
    Rgrad_likeli_L = zeros([length(thetaL) save_it]);
    Rgrad_const_L = zeros([length(thetaL) save_it]);
    
    Rdelta_L = zeros([length(thetaL) save_it]);
    
    Rconst = zeros([1 save_it]);
end;

Z=sum(exp(thetaL(gen_len) + sum(theta1(gen_subs),2)))/gen_samp; %normalization

grad_likeli_aa=data_aa_f; %the aa part of the log likelihood gradient -fixed
grad_likeli_L=data_dlen;

%%
delta=Inf;
i=0;
while delta>tol && i<max_it %GD loop continues until selection factor change is small enough
    i=i+1;
    %%
    theta1_c=theta1;
    thetaL_c=thetaL;
        
        %multiplying the correct qs for each sequence
        Q_gen=Z^-1*exp(thetaL(gen_len) + sum(theta1(gen_subs),2));
        %sum over the products of all sequences with same aa and position combination
        grad_const_aa=accumarray(gen_inds,Q_gen(gen_numseq),[max_aa,maxL,maxL])/gen_samp + pc;       
        grad_const_L=accumarray(gen_len,Q_gen,[maxL 1])/gen_samp;
        %%
        norm_fac = 1; %norm_fac=norm([grad_const_L(:) ; grad_const_aa(:)]); %norm of the entire constrain gradient "vector"
        proj = 1; %proj=sum(grad_likeli_L(:).*grad_const_L(:))/norm_fac+sum(grad_likeli_aa(:).*grad_const_aa(:))/norm_fac; %projection of the likelihood gradient in the direction of the constrain gradient (normalized)
        
        delta_L=GD_step_L*(grad_likeli_L-proj*grad_const_L/norm_fac);
        delta_aa=GD_step_aa*(grad_likeli_aa-proj*grad_const_aa/norm_fac); %gradient assent change that keeps the constrain
        delta_aa(:,:,~Lrange)=0; %set all changes in low statistics zone to zero
                
        theta1=theta1+delta_aa; %update qs
        thetaL=thetaL+delta_L;
    %disp(['finished max iteration ',num2str(i) ,', updating expectation']);

    D1=theta1(isfinite(theta1))-theta1_c(isfinite(theta1_c));
    DL=thetaL(isfinite(thetaL))-thetaL_c(isfinite(thetaL));
    delta=sqrt(nansum(DL.^2)+sum(D1.^2)); 
    %disp(['delta=' num2str(delta)]);
    
    if mod(i,M)==0
        disp(['max iteration # ' num2str(i) ' **** delta=' num2str(delta)]); 
        if his
            %         saving history of the run
            ii=i/M;
            Rtheta1(:,:,:,ii)=theta1;
            RthetaL(:,ii)=thetaL;
            Rgrad_likeli_aa(:,:,:,ii)=grad_likeli_aa;
            Rgrad_const_aa(:,:,:,ii)=grad_const_aa;
            Rgrad_likeli_L(:,ii)=grad_likeli_L;
            Rgrad_const_L(:,ii)=grad_const_L;
            Rdelta_aa(:,:,:,ii)=delta_aa;
            Rdelta_L(:,ii)=delta_L;
            Rconst(ii)=mean(Q_gen);
        end;
    end;
    
end;

theta1=bsxfun(@minus,theta1,log(sum(exp(theta1).*bsxfun(@rdivide,gen_aa_f,shiftdim(gen_dlen,-2))))); %guage fixing
theta1(gen_aa_f==0)=0; %combinations that don't appear at all in gen have selection factors exactly 1 (for lack of data)
qL=exp(thetaL);
Z=sum(qL(gen_len).*exp(sum(theta1(gen_subs),2)))/gen_samp;

if his
    save(output_name,'theta1','Z','i','qL','R*');
else
    save(output_name,'theta1','Z','i','qL');
end;
disp('done');
