function [xnew,MemoryNode_new,MemoryCluster_new,powerUO,LinkMat,varargout]= ...
EsaR_iter(tri_actif,nontri,xinit,y,V,MemoryNode,MemoryCluster,...
Dist,LinkMatOld,alpha, proba_indep, type_fail,varargin)
%
% usage [xnew,MemoryNode_new,MemoryCluster_new,powerUO,LinkMat,varargout]=...
% EsaR_iter(tri_actif,nontri,xinit,y,V,MemoryNode,MemoryCluster,...
% Dist,LinkMatOld,alpha,etat, proba_indep, type_fail,varargin)
%
% Input : tri_actif : ntrix3 matrix of nodes that are in the embedded triangles
% nontri : 1:n \ tri_actif : nodes that are in the network but
% not in the embedded trianlges
% xinit : current posterior mean
% y : observations/sigma (normalised observations)
% V : matrix to be inverted
% MemoryNode : memory that each node has of its neighbors. This
% is used in the robust update
% MemoryCluster : inside an embedded triangle, each sensor i keeps
% in memory the value \hat{y_k}^m for k = the two
% other nodes that in the same triangle as i
%
% Dist : matrix with Distance between nodes in the network,
% when these nodes are linked
% LinkMatOld : Matrix giving the state of each link.
% LinkMatOld(i,j) = 1 if link is active, =1 if it is
% the first iteration where link is inactive =-2 if
% it is the second step, etc,...
% LinkMatOld(i,j) = 0 if there are no link between i
% and j
% alpha : exponent in computation of the Energy as a fct of Dist
% proba_indep = [p1 M], with p1 probability of link failure at each iteration
% and M maximum number of iteratinons where a link or
% node can be off.
% type_fail = 'node_fail' or 'link_fail'
% varargin{1} = nodeOff = indices of nodes that are sleeping (for
% type_fail = node_fail)
%
% Output : xnew : posterior mean at next iteration
% MemoryNode_new: updated MemoryNode
% MemoryCluster_new : updated MemoryCluster
% powerUO : power unidirectional and omni consumed at this
% iteration
% LinkMat : new matrix indicating links that are ON or OFF (and
% if OFF, at which stage)
% varargout = updated nodeOFF : indices of the new nodes that
% are sleeping
%
% called by Esa_robust.m
% Adjacency matrix
dag = (LinkMatOld ~= 0);
n=size(V,1);
taille = size(tri_actif);
nbtri=taille(1); % number of disjoint embedded triangles
nbnodes = taille(1)*taille(2);
DD=spdiags(sparse(V),0);
J=sparse(diag(DD));
for k=1:nbtri
trik = tri_actif(k,:);
J(trik,trik) = V(trik,trik);
end;
K = J-V;
x=xinit;
MemoryCluster_new = MemoryCluster;
elemtri =sort(reshape(tri_actif, taille(1)*taille(2),1));
switch type_fail
case 'link_fail'
[LinkMat] = matfailureN_L(LinkMatOld, proba_indep,type_fail);
case 'node_fail'
nodeOff=varargin{1};
[LinkMat, ind_fail] = matfailureN_L(LinkMatOld, proba_indep,...
type_fail,nodeOff);
varargout(1)={ind_fail};
otherwise
error('type_fail in EsaR_iter should be link_fail or node_fail')
end
% product Ktri_x = Ktri*x but now with link failure and/or node down
% new value transmitted :
% NOTE : linkmat(i,j) gives the transmission node from i to j
% here what we need is transmission from x_atnode to x_trans
% hence we need to transpose _linkmat_
% from LinkMat build a matrix of 0 and 1 :
p1 = proba_indep(1);
linkmat = (LinkMat> 0);
% x_trans = x transmitted : values that were transmitted
x_trans = (K.*linkmat')*x;
Kfailed = K.*(1 - linkmat');
% x_nontrans = x non transmitted : values that were NOT transmitted
x_nontrans = sum(Kfailed.*MemoryNode,2);
Kx = x_trans + x_nontrans;
for k=1:nbtri % For all embedded triangles
trik=tri_actif(k,:);
Jtrik=J(trik,trik);
% UPDATE STEP : either ynewk is computed from x_transmitted and
% x_nontransmitted.
ynewk = y(trik) + Kx(trik);
for indtri = 1:3
memk = MemoryCluster{indtri}(k,:);
one = indtri;
switch type_fail
case 'link_fail'
oneOK =1;
otherwise
if length(ind_fail)>0
oneOK = ~sum(ind_fail == trik(one));
else
oneOK=1;
end
end % end switch type_fail
switch oneOK
case 0 % node 'one' not solved
x(trik(one)) = xinit(trik(one));
case 1
twothreeind = setdiff(1:3,indtri);
two= (twothreeind(1));
three = (twothreeind(2));
% If the link between two and one or between three and one
% fails, we use the value in the memory
if ~linkmat(trik(two), trik(one))
ynewk(two) = memk(1);
end
if ~linkmat(trik(three), trik(one))
ynewk(three) = memk(2);
end
% SOLVE step when transmission was successful
MemoryCluster_new{indtri}(k,:) = [ynewk(two),ynewk(three)];
xtmp = Jtrik\ynewk;
x(trik(one))= xtmp(one);
end% end switch oneOK
end % end for indtri
end % for k
% If a node if off, it will not receive its new value, even if
% its neighbors try to communicate with him
if strcmp(type_fail,'node_fail')
% node who fails do not update their value
if length(ind_fail) >0
x(ind_fail) = xinit(ind_fail);
end
end
% Update for other nodes:
% Gauss-seidel type: we use the newly computed values for x
% Ktri_x = Ktri*x;
x_trans = (K.*linkmat')*x;
Kfailed = K.*(ones(size(linkmat)) - linkmat');
x_nontrans = sum(Kfailed.*MemoryNode,2);
Kx = x_trans + x_nontrans;
D = spdiags(V(nontri,nontri),0);
if length(nontri) >0
x(nontri,1) = (y(nontri,1) + Kx(nontri,1))./D;
end
% Again, if some node in nontri are OFF, their value is not updated
%---- Now this part is treated above ...
if ~strcmp(type_fail,'link_fail')
% node who are sleeping do not update their value
if length(ind_fail) >0
x(ind_fail) = xinit(ind_fail);
end
end
xnew = x;
% Total energy spent by the sensors
% Here you take into account if the sleeping modes of the sensors
% are synchronised or not synchronised
if strcmp(type_fail,'node_fail')
if length(ind_fail) > 0
powerUO = TrianPowerHeadR(tri_actif,Dist,J,K,alpha,ind_fail);
else
powerUO = TrianPowerHead(tri_actif,Dist,J,K,alpha);
end
else % type_fail = 'link_fail'
powerUO = TrianPowerHead(tri_actif,Dist,J,K,alpha);
end
%x_atnode is the value transmitted at this iteration: hence
% if the transmission has been successful it must be recorded
MemoryNode_new = newMemory(linkmat, dag, x, MemoryNode);
%------------------------------------------------------------------