[EC,ER]=CRcheq(C,R,P)

function [EC,ER]=CRcheq(C,R,P)
% function [EC,ER]=CRcheq(C,R,P)
% check if C, R are solutions of the (adversarial or drunk) optimality equations
%
% C:	cop cost matrix
% R:	robber cost matrix (use empty matrix for drunk optim. equations)
% E:	satisfaction matrix
%		EC(mc,mr)=0 iff (mc,mr) cop eq. is satisfied
%		ER(mc,mr)=0 iff (mc,mr) rob eq. is satisfied
NC=size(C);
NR=size(R);

N=NC(1);
if	max(NR)>0  & length(NC)==2			% adv.  robber, one cop
	[EC,ER]=CRcheqA2(C,R,P,N);
elseif	max(NR)==0 & length(NC)==2			% drunk robber, one cop
	[EC,ER]=CRcheqD2(C,R,P,N);
elseif	max(NR)>0  & length(NC)==3			% adv.  robber, two cops
	[EC,ER]=CRcheqA3(C,R,P,N);
elseif	max(NR)==0 & length(NC)==3			% drunk robber, two cops
	[EC,ER]=CRcheqD3(C,R,P,N);
else
	EC=[];
	ER=[];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [EC,ER]=CRcheqA2(C,R,P,N)
EC=zeros(N,N);
ER=zeros(N,N);
PP=P+eye(N);
for mc=1:N
	for mr=1:N
		if 	mc==mr
			if C(mc,mr)==0; EC(mc,mr)=0; end
			if R(mc,mr)==0; ER(mc,mr)=0; end
		else
			Nc=find(PP(mc,:)>0);
			minR=min(R(Nc,mr));
			EC(mc,mr)=abs(C(mc,mr)-1-minR);
			Nr=find(PP(mr,:)>0);
			maxC=max(C(mc,Nr));
			ER(mc,mr)=abs(R(mc,mr)-maxC);
		end
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [EC,ER]=CRcheqD2(C,R,P,N)
EC=zeros(N,N);
ER=[];
PP=P+eye(N);
for mc=1:N
	for mr=1:N
		if 	mc==mr
			if C(mc,mr)==0; EC(mc,mr)=0; end
		else
			Nc=find(PP(mc,:)>0);
			CC=ones(1,N)*1e15;
			for nc=Nc
				CC(nc)=1;
				P1=P; P1(nc,:)=0; P1(:,nc)=0;
				Nr=find(P1(mr,:)>0);
				for nr=Nr
					CC(nc)=CC(nc)+P1(mr,nr)*C(nc,nr);
				end
			end
			minC=min(CC);
			EC(mc,mr)=abs(C(mc,mr)-minC);
		end
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [EC,ER]=CRcheqA3(C,R,P,N)
EC=zeros(N,N,N);
ER=zeros(N,N,N);
PP=P+eye(N);
for mc1=1:N
	for mc2=1:N
		for mr=1:N
			if 	mc1==mr || mc2==mr
				if C(mc1,mc2,mr)==0; EC(mc1,mc2,mr)=0; end
				if R(mc1,mc2,mr)==0; ER(mc1,mc2,mr)=0; end
			else
				Nc1=find(PP(mc1,:)>0);
				Nc2=find(PP(mc2,:)>0);
				minR=min(min(R(Nc1,Nc2,mr)));
				EC(mc1,mc2,mr)=abs(C(mc1,mc2,mr)-1-minR);
				Nr=find(PP(mr,:)>0);
				maxC=max(C(mc1,mc2,Nr));
				ER(mc1,mc2,mr)=abs(R(mc1,mc2,mr)-maxC);
			end
		end
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [EC,ER]=CRcheqD3(C,R,P,N)
EC=zeros(N,N,N);
ER=[];
PP=P+eye(N);
for mc1=1:N
	for mc2=1:N
		for mr=1:N
			if 	mc1==mr || mc2==mr
				if C(mc1,mc2,mr)==0; EC(mc1,mc2,mr)=0; end
			else
				Nc1=find(PP(mc1,:)>0);
				Nc2=find(PP(mc2,:)>0);
				CC=ones(N,N)*1e15;
				for nc1=Nc1
					for nc2=Nc2
						CC(nc1,nc2)=1;
						P1=P; P1(nc1,:)=0; P1(:,nc1)=0; P1(nc2,:)=0; P1(:,nc2)=0;
						Nr=find(P1(mr,:)>0);
						for nr=Nr
							CC(nc1,nc2)=CC(nc1,nc2)+P1(mr,nr)*C(nc1,nc2,nr);
						end
					end
				end
				minC=min(min(CC));
				EC(mc1,mc2,mr)=abs(C(mc1,mc2,mr)-minC);
			end
		end
	end
end