
function [ynew,P,R,prf]=learning_BIS(std_P,mean_P,y,R,P,Mshall,epsilon,A1,A2,A3,neq,nsh,truenonzero,truezeros,njump,gain,picknonzero,picknojumps);
% learning based on solution 0=A1*Pi2*PHI2^2+A2*Pi2*PHI2+A3*Pi2;
% Then 0=A1*Pi2*PHI2*y(t)+A2*y(t)+A3*y(t-1);
% where PLM=> E y(t+1)=Pi2*PHI2*y(t)

R0=R;
P0=P;
prf=0;


%% ALM
% note that the order of variables has changed, Thus the reordering of
% Mshall relative to RE solution(end:-1:1,end:-1:1)
ynew=(A1*P+A2)\(-A3*y+Mshall*epsilon);
%% update parameters
Pstates=P(njump,truenonzero);
X=y(truenonzero);% states
% X(nsh+1:end)=X(nsh+1:end)+randn(length(X)-nsh)/100;% measurement error
R=R+gain*(X*X'-R); %timing!!
% forcmx=kron_d_f(eye(length(njump)),R\X);
errort=ynew(njump)-Pstates*X;%

Pstates=Pstates+gain*(R\(X*errort'))';
% PP=vec((Pstates)')+gain*forcmx*errort;
% Pstates=reshape(PP,size((Pstates)'))';
P=picknojumps*Pstates*picknonzero;

%% Projection Facility: including non invertibility (of ALM, Covariance of regressors and unusual estimates (x*standard errors)

if rcond(R)<1e-5
        R=R0+0.00001*eye(size(R0,1));
%     P=P0;
%     0.99*diag(rand(size(P0,1),1));
    prf=1;
    return
end
MMM=(abs(mean_P));

% if any(abs((ynew(truenonzero))))>10*(diag(R).^0.5)
try
if any(ynew(truenonzero)>MMM+20*((diag(R)).^0.5))
    P=P0;
    prf=1;
    R=R0+0.00001*eye(size(R0,1));
end
catch
       P=P0;
    prf=1;
    R=R0+0.00001*eye(size(R0,1));
end
try
if max(abs(eig((A1*P+A2)\A3)))>1
        R=R0;
        P=P0;
        %0.99*diag(rand(size(P0,1),1));
        prf=1;
        return
        
end
catch
      R=R0;
        P=P0;
        %0.99*diag(rand(size(P0,1),1));
        prf=1;
end


% if rcond((A1*P+A2))<1e-8||rcond(R)<1e-8||any(abs(vec(Pstates))>20*abs(std_P));
%     R=R0;
%     P=mean_P;
%     prf=1;
%     return
% end
% try % sometimes it occurred that the previous was passed and the inverse had NAN or INF
%     if max(abs(eig((A1*P+A2)\A3)))>1
%         R=R0;
%         P=mean_P;%+0.01*diag(rand(size(P0,1),1));
%         prf=1;
%         return
%     end
% catch
%     R=R0;
%     P=mean_P;%+0.01*diag(rand(size(P0,1),1));
%     prf=1;
%     return
% end
% tmp2=eye(nsh);
% Pi2P=[tmp2(:,end:-1:1),zeros(nsh,neq-nsh);
%     P(nsh+1:neq,:)];


%
% st=P(truenonzero,truenonzero)*st+Mshall(truenonzero,:)*epsilon;

% Ainv=tmp\eye(neq);




end