function W = solver(B)
qpz8U = 1000;
rand('seed',32);
mfh8f = [];
W = zeros(0,4);
[joN2B,S_Hl6] = meshgrid(1:size(B,2),1:size(B,1));
qGeaN = size(B,1);
BrRuc = size(B,2);
Vx8cb = zeros(size(B));
OpRdC = qGeaN*BrRuc;
Vx8cb(1:OpRdC) = 1:OpRdC;
ZDXni =  [ 1000*ones(qGeaN,1) , Vx8cb(:,1:end-1)];
kx3GI =  [ Vx8cb(:,2:end) , 1000*ones(qGeaN,1) ];
ZcOYj =  [ Vx8cb(2:end,:) ; 1000*ones(1,BrRuc) ];
rPWRX =  [ 1000*ones(1,BrRuc) ; Vx8cb(1:end-1,:)];
dov7h = zeros(qGeaN*BrRuc,4);
dov7h(:,1:2) = 1;
[rM7r8,l0h8r] = sort([ ZDXni(:) , kx3GI(:) , ZcOYj(:),rPWRX(:) ],2);
dov7h = dov7h(sub2ind([qGeaN*BrRuc,4],repmat((1:qGeaN*BrRuc)',1,4),l0h8r));
[ZDXni,kx3GI] = ind2sub([qGeaN*BrRuc,4],find(rM7r8==1000) );
rM7r8(sub2ind([qGeaN*BrRuc,4],ZDXni,kx3GI)) = ZDXni;
dNgvO = repmat((1:qGeaN*BrRuc)',1,4);
D_R4R = 1;
fED4B = find(B);
ubppn = sortrows([B(fED4B(:)),fED4B]);
fED4B = ubppn(:,2);
[Values,khaoF,JL_q5] = unique(ubppn(:,1),'first');
eqsIY = [khaoF(2:end)-1 ; numel(fED4B)];
mfJBd = histc(ubppn(:,1),Values);
Yafxv = numel(Values);
Y2csy =  - ( mfJBd .^ 1 ) .* ( Values.^ 4 );
emMtU = (mfJBd .* Values);
MhJJY = qGeaN*BrRuc;
slE3R = B>0;
n3_2T = slE3R;
X8JIt = slE3R;
k_hOA = zeros(Yafxv,1);
W = [];
while 1
zyY04 = n3_2T;
Zl3AU = X8JIt;
[F_kiP,AtOVQ] = fBute(slE3R);
if isempty(F_kiP)
break
end
if isempty(find(AtOVQ(:,5)>0,1))
break
end
[Gwod_,l0h8r] = max(AtOVQ(:,5));
k_hOA(AtOVQ(l0h8r,6)) = 1;
fK0vE = F_kiP(AtOVQ(l0h8r,1):AtOVQ(l0h8r,2),:);
W = [W ; fK0vE];
[ slE3R , X8JIt , n3_2T ] = ...
JdjaE( slE3R , X8JIt , n3_2T , fK0vE );
end
if isempty(W)
W=zeros(0,4);
end
function [tUR7l,rZt_Y,FztNu] = JdjaE(slE3R,sNcLC,DACpK,eSFgN)
tUR7l = slE3R;
rZt_Y = sNcLC;
FztNu = DACpK;
kXBgc = eSFgN(:,1)==eSFgN(:,3) & eSFgN(:,2)==eSFgN(:,4);
l0h8r = sub2ind( size(slE3R) , eSFgN(:,[1 3]) , eSFgN(:,[2 4]) );
tUR7l(l0h8r) = 1;
TRQWU = eSFgN(:,1) ~= eSFgN(:,3);
rZt_Y(l0h8r(TRQWU,:)) = 1;
TRQWU = eSFgN(:,2) ~= eSFgN(:,4);
FztNu(l0h8r(TRQWU,:)) = 1;
end
function [ExWnQ,AWXRW] = fBute(lcO8r)
AWXRW = [];
ExWnQ = [];
[Gwod_,QM0nt] = sort(Y2csy);
VsKIM = sum(emMtU(k_hOA==0));
NlwgR = VsKIM/sum(emMtU);
fcAey = numel(find(lcO8r==0));
D_R4R = VsKIM / fcAey * .3;
for BDXiZ=1:Yafxv
fCh9l = QM0nt(BDXiZ);
if khaoF(fCh9l)<eqsIY(fCh9l) && ~k_hOA(fCh9l);
wbxjY = size(ExWnQ,1);
[ddcek,SxmsU]=GPgeP(lcO8r,fCh9l);
if isempty(SxmsU)
continue
end
SxmsU(:,1:2) = SxmsU(:,1:2) + wbxjY;
ExWnQ = [ExWnQ ; ddcek];
AWXRW = [AWXRW ; [SxmsU , ones(size(SxmsU,1),1)*fCh9l ] ];
if ~isempty(ExWnQ)
break
end
end
end
end
function [ExWnQ,AWXRW] = oMXBu(lcO8r)
AWXRW = [];
ExWnQ = [];
for fCh9l=1:Yafxv
if khaoF(fCh9l)<eqsIY(fCh9l) && ~k_hOA(fCh9l);
wbxjY = size(ExWnQ,1);
[ddcek,SxmsU]=GPgeP(lcO8r,fCh9l);
if isempty(SxmsU)
continue
end
SxmsU(:,1:2) = SxmsU(:,1:2) + wbxjY;
ExWnQ = [ExWnQ ; ddcek];
AWXRW = [AWXRW ; [SxmsU , ones(size(SxmsU,1),1)*fCh9l ] ];
end
end
end
function [eJoid,i6p9H] = GPgeP(lcO8r,cfZeU)
i6p9H = [];
YrzhL = khaoF(cfZeU);
jEFr2 = eqsIY(cfZeU);
GgR44 = jEFr2-YrzhL+1;
Jb69O = fED4B(YrzhL:jEFr2);
mEfnu = 1001*ones(qGeaN*BrRuc,GgR44);
DmDJt = Jb69O;
uLTXj = (1:GgR44)';
HIprO = sub2ind([qGeaN*BrRuc,GgR44],DmDJt,uLTXj);
W63zV = min( qGeaN+BrRuc , ceil(Values(cfZeU)*GgR44/2) );
mEfnu(HIprO) = 0;
while ~isempty(HIprO) && ~isempty(find(mEfnu(HIprO) < W63zV ,1,'first'))
iolNY = rM7r8(DmDJt,:);
U09Yc = B(iolNY)~=Values(cfZeU);
p56V4 = ones(size(iolNY)) + 1000 * ...
( U09Yc & ( B(iolNY)>0 | ...
(zyY04(iolNY) & dov7h(DmDJt,:)) ...
| (Zl3AU(iolNY) & ~dov7h(DmDJt,:) ) ) ) + ...
qpz8U * (lcO8r(iolNY) & U09Yc);
aYW21 = p56V4 + repmat(mEfnu(HIprO),1,4);
[ZDXni,kx3GI] = ind2sub( size(aYW21), (1:prod(size(aYW21)))' );
vxizI = sub2ind([qGeaN*BrRuc,GgR44] , iolNY , repmat(uLTXj,1,4));
ubppn = sortrows([ vxizI(:) , aYW21(:)],[1 2]);
[aqmKz , SXRde] = unique(ubppn(:,1),'first');
OqE5H = ubppn(SXRde,2);
CUzwa = OqE5H<1000 & OqE5H<mEfnu(aqmKz);
HIprO = aqmKz(CUzwa);
mEfnu(HIprO) = min( mEfnu(HIprO) , OqE5H(CUzwa));
[DmDJt,uLTXj] = ind2sub([qGeaN*BrRuc,GgR44] , HIprO);
end
[j3jhe,sCrXN] = sort(cumsum(sort(mEfnu,2),2),1);
obIE3 = min(4,size(j3jhe,1));
ARgLC = repmat( (2:GgR44)*Values(cfZeU) , obIE3 ,1);
WaqBq = repmat( (2:GgR44)*Values(cfZeU) , obIE3 ,1) - ...
j3jhe(1:obIE3,2:end) ./ ( min(1,1e3*ARgLC.^(1)) ) * (1+D_R4R);
[gYyYF,tFeOq] = sort(WaqBq(:),'descend');
Oe3DO = min(3,numel(tFeOq));
tFeOq = tFeOq(1:Oe3DO);
gYyYF = gYyYF(1:Oe3DO);
tFeOq = tFeOq(gYyYF(:)>0);
[UeL50,lAi6p] = ind2sub([obIE3 , GgR44-1],tFeOq);
COcyB = sCrXN(sub2ind( size(sCrXN) , UeL50 , lAi6p+1));
uX3ig = lAi6p + 1;
eJoid = [];
if isempty(COcyB)
return
end
for nDcOz = 1:numel(COcyB)
fhm4L = B7dl6(COcyB(nDcOz),uX3ig(nDcOz));
Nu_3R(nDcOz) = size(eJoid,1) + 1;
eJoid = [eJoid ; fhm4L];
iVPez(nDcOz) = size(eJoid,1);
hEQki(nDcOz) = sum(fhm4L(:,1)==fhm4L(:,3) & fhm4L(:,2)==fhm4L(:,4));
bATew(nDcOz) = size(fhm4L,1) - hEQki(nDcOz);
ET5DX = bATew(nDcOz) + qpz8U * hEQki(nDcOz);
qB_mj = uX3ig(nDcOz) * Values(cfZeU);
MEfYu(nDcOz) = qB_mj - ET5DX;
end
[Gwod_,l0h8r] = max(MEfYu);
i6p9H = [Nu_3R(:),iVPez(:),hEQki(:),bATew(:),MEfYu(:)];
function PuyDN = B7dl6(jk3Nc,mWWYL)
[Gwod_,FcBXw] = sort(mEfnu(jk3Nc,:)');
FcBXw = FcBXw(1:mWWYL);
HIprO = jk3Nc;
cWo7H = FcBXw;
eQ4GT = ones(mWWYL,1);
KvoW2 = HIprO(eQ4GT);
eYvwr = zeros(0,1);
Etjni = zeros(0,1);
nMvMk = zeros(0,1);
while ~isempty(cWo7H)
XMi9V = mEfnu(sub2ind(size(mEfnu),KvoW2(:),cWo7H));
XoJ_i = rM7r8(KvoW2,:);
DYu7t = mEfnu( sub2ind(size(mEfnu),XoJ_i , repmat(cWo7H(:),1,4) ) );
XBjAN = repmat(XMi9V,1,4);
FyMhL = XBjAN - DYu7t == (qpz8U+1);
eha1S = DYu7t<XBjAN;
[jan3a,Gwod_,L95b1] = unique(XoJ_i(:));
Count = histc(XoJ_i(eha1S),jan3a);
[ngo1v , gFhI3 ] = max(eha1S .* reshape(Count(L95b1),size(eha1S)) , [] , 2);
meT2l = ngo1v>0;
rPybY = find(meT2l);
T3GlN = sub2ind( size(eha1S) , rPybY , gFhI3(meT2l) );
Bv4oD = FyMhL(T3GlN);
cWo7H = cWo7H(rPybY);
Yo1fU = KvoW2(rPybY);
eYvwr = [eYvwr ; Yo1fU];
KvoW2 = XoJ_i(T3GlN);
Etjni = [Etjni ; KvoW2];
EFZFn = Yo1fU(Bv4oD);
if ~isempty(EFZFn)
nMvMk = [nMvMk ; EFZFn];
end
end
pin0h = min(eYvwr,Etjni);
pAnIJ = max(eYvwr,Etjni);
[eYvwr,Etjni] = ...
ind2sub([qGeaN*BrRuc,qGeaN*BrRuc], unique(sub2ind([qGeaN*BrRuc,qGeaN*BrRuc],pin0h,pAnIJ)) );
nMvMk = unique(nMvMk);
[ZDXni,kx3GI] = ind2sub([qGeaN,BrRuc],eYvwr);
[ZcOYj,rPWRX] = ind2sub([qGeaN,BrRuc],Etjni);
[j6MTz,MWueN] = ind2sub([qGeaN,BrRuc],nMvMk);
PuyDN = [ZDXni,kx3GI,ZcOYj,rPWRX];
PuyDN = [PuyDN; [j6MTz,MWueN,j6MTz,MWueN]];
end
end
end