0001 function [F_rfc,F_rfc_z,T] = smctp2arfm(P,F,c_m,SideInfo)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069 T(1,:)=clock;
0070
0071
0072
0073 ni = nargin;
0074 no = nargout;
0075 error(nargchk(2,4,ni));
0076
0077 if ni < 3
0078 c_m=[];
0079 end
0080
0081 if ni < 4
0082 SideInfo = [];
0083 end
0084
0085 if isempty(c_m)
0086 c_m=1;
0087 end
0088 if isempty(SideInfo)
0089 SideInfo=0;
0090 end
0091
0092
0093
0094 Zstr = '123456789'; Zstr=Zstr';
0095
0096 r = length(P);
0097 n = length(F{1,1});
0098
0099
0100
0101 P = mat2tmat(P);
0102
0103 T(2,:)=clock;
0104
0105
0106
0107
0108 for i = 1:r
0109 QQ{i,1} = F{i,1};
0110 QQ{i,1} = mat2tmat(QQ{i,1},1);
0111 end
0112
0113 T(3,:)=clock;
0114
0115
0116
0117
0118 for i = 1:r
0119
0120 if isempty(F{i,2})
0121 QQ{i,2} = F{i,1};
0122 QQ{i,2} = QQ{i,2}';
0123 else
0124 QQ{i,2} = F{i,2};
0125 end
0126
0127 QQ{i,2} = mat2tmat(QQ{i,2},-1);
0128
0129 end
0130
0131 T(4,:)=clock;
0132
0133
0134
0135 Q = zeros(n*r,n*r);
0136 I = 0:r:(n-1)*r;
0137 for z=1:r
0138 Q0 = kron(QQ{z,1},P);
0139 Q(I+z,:) = Q0(I+z,:);
0140 end
0141
0142 T(5,:)=clock;
0143
0144
0145
0146 Qh = zeros(n*r,n*r);
0147 I = 0:r:(n-1)*r;
0148 for z=1:r
0149 Q0 = kron(QQ{z,2},P);
0150 Qh(I+z,:) = Q0(I+z,:);
0151 end
0152
0153 T(6,:)=clock;
0154
0155
0156
0157
0158 Qt = Q*Qh;
0159 ro = mc2stat(Qt(1:r*(n-1),1:r*(n-1)));
0160 ro = [ro zeros(1,r)];
0161
0162
0163
0164
0165 Qth = Qh*Q;
0166 roh = mc2stat(Qth(r+1:r*(n),r+1:r*(n)));
0167 roh = [zeros(1,r) roh];
0168
0169 T(7,:)=clock;
0170
0171
0172
0173
0174 FF = Q.*(ro'*ones(1,n*r)) + Qh.*(roh'*ones(1,n*r));
0175
0176
0177
0178
0179 I1 = find(ro>0); I2 = find(ro<=0);
0180 ro_inv = ro; ro_inv(I1) = 1./ro(I1); ro_inv(I2) = zeros(1,length(I2));
0181 Qr = Qh' .* (ro_inv'*roh);
0182
0183
0184 I1 = find(roh>0); I2 = find(roh<=0);
0185 roh_inv = roh; roh_inv(I1) = 1./roh(I1); roh_inv(I2) = zeros(1,length(I2));
0186 Qrh = Q' .* (roh_inv'*ro);
0187
0188
0189
0190
0191 FF1 = Qr.*(ro'*ones(1,n*r)) + Qrh.*(roh'*ones(1,n*r));
0192
0193 T(8,:)=clock;
0194
0195
0196
0197 F_rfc = zeros(n,n);
0198 EYE = eye(n*r,n*r); Eye1 = eye(r,r);
0199 Zero1 = zeros(r,1); Zero2 = zeros(r,r);
0200 One1 = ones(r,1); One2 = ones(r,r);
0201
0202 if SideInfo == 1
0203 [F_rfc_z{1:r,1:r}] = deal(zeros(n,n));
0204 end
0205
0206
0207
0208
0209
0210 for i=1:n-1
0211
0212
0213 for j=i+1:n
0214
0215 I = r*(i-1)+1:r*i;
0216 J = r*(j-1)+1:r*j;
0217 I1 = r*(i+1-1)+1:r*(i+1);
0218 J1 = r*(j-1-1)+1:r*(j-1);
0219 I1J1 = r*(i+1-1)+1:r*(j-1);
0220
0221 x = FF(I,J);
0222 y = FF(J,I);
0223
0224 Min = sum(FF(I,:),2)';
0225 Max = sum(FF(J,:),2)';
0226
0227 Ro = ro(I);
0228 Roh = roh(J);
0229
0230
0231 Min = Ro; Max = Roh;
0232
0233 Min = Min'*One1';
0234 Max = Max'*One1';
0235
0236
0237
0238 C = y'./Min; Ch = x'./Max;
0239 E = 1-sum(y'./Min,2); Eh = 1-sum(x'./Max,2);
0240
0241
0242 if j-i >= 2
0243
0244
0245
0246 d1 = Qr(I,I1)*(1-sum(Qrh(I1,I),2));
0247 D = d1;
0248 d1h = Qrh(J,J1)*(1-sum(Qr(J1,J),2));
0249 Dh = d1h;
0250 d0 = sum(Qr(I,I1J1),2);
0251 E = E-d0;
0252 d0h = sum(Qrh(J,I1J1),2);
0253 Eh = Eh-d0h;
0254 r1 = Qr(I,I1)*Qrh(I1,I);
0255 R = r1;
0256 r1h = Qrh(J,J1)*Qr(J1,J);
0257 Rh = r1h;
0258
0259
0260
0261 d3 = sum(Q(I,I1J1),2);
0262 E3 = 1-d3;
0263 d3h = sum(Qh(J,I1J1),2);
0264 E3h = 1-d3h;
0265
0266 else
0267
0268 D = Zero1; Dh = Zero1;
0269 R = Zero2; Rh = Zero2;
0270 E3 = One1; E3h = One1;
0271
0272 end
0273
0274
0275 if j-i >= 3
0276
0277 Eye = EYE(1:r*(j-i-2),1:r*(j-i-2));
0278
0279
0280
0281 II = r*(i+1-1)+1:r*(j-2);
0282 JJ = r*(i+2-1)+1:r*(j-1);
0283
0284 A = Qr(II,JJ);
0285 Ah = Qrh(JJ,II);
0286 Inv = inv(Eye-A*Ah);
0287
0288 a = Qr(I,JJ);
0289 ah = Qrh(J,II);
0290 b = Qr(II,J);
0291 bh = Qrh(JJ,I);
0292
0293 e = 1 - sum(Qr(II,r*(i+2-1)+1:r*(j)),2);
0294 eh = 1 - sum(Qrh(JJ,r*(i-1)+1:r*(j-2)),2);
0295
0296 C = C + a*Ah*Inv*b;
0297 Ch = Ch + ah*Inv*A*bh;
0298 D = D + a*eh+a*Ah*Inv*A*eh;
0299 Dh = Dh + ah*Inv*e;
0300 E = E + a*Ah*Inv*e;
0301 Eh = Eh + ah*Inv*A*eh;
0302 R = R + a*bh+a*Ah*Inv*A*bh;
0303 Rh = Rh + ah*Inv*b;
0304
0305
0306
0307 A3 = Q(II,JJ);
0308 A3h = Qh(JJ,II);
0309 Inv3 = inv(Eye-A3*A3h);
0310 c3 = Q(I,JJ);
0311 c3h = Qh(J,II);
0312 e3 = 1 - sum(Q(II,r*(i+2-1)+1:r*(j-1)),2);
0313 e3h = 1 - sum(Qh(JJ,r*(i+1-1)+1:r*(j-2)),2);
0314
0315 E3 = E3 + c3*A3h*Inv3*e3;
0316 E3h = E3h + c3h*Inv3*A3*e3h;
0317
0318 end
0319
0320 if ~(i == 1 & j == n)
0321
0322
0323 if j == n
0324 pS1 = Zero1; pS2=Zero2; pS3=Zero1;
0325 else
0326
0327
0328 AS = (Eye1-Rh)\Dh;
0329 ES = E + C*AS;
0330 BS = (Eye1-Rh)\Ch;
0331 RS = R + C*BS;
0332 pS1 = (Eye1-RS)\ES;
0333
0334
0335 pS2 = (Eye1-Rh)\Ch;
0336
0337
0338 pS3 = E3h;
0339
0340 end
0341
0342
0343 if i == 1
0344 pH1 = Zero1; pH2=Zero2; pH3=Zero1;
0345 else
0346
0347
0348 AH = (Eye1-R)\D;
0349 EH = Eh + Ch*AH;
0350 BH = (Eye1-R)\C;
0351 RH = Rh + Ch*BH;
0352 pH1 = (Eye1-RH)\EH;
0353
0354
0355 pH2 = (Eye1-R)\C;
0356
0357
0358 pH3 = E3;
0359 end
0360
0361 else
0362
0363
0364 pS1 = One1;
0365 pS2 = Eye1;
0366 pS3 = E3h*sum(Roh)/sum(Ro+Roh);
0367
0368
0369 pH1 = One1;
0370 pH2 = Eye1;
0371 pH3 = E3*sum(Ro)/sum(Ro+Roh);
0372
0373 end
0374
0375
0376 F_rfc(i,j) = Roh*diag(pS3)*pS2*pS1;
0377
0378
0379 F_rfc(j,i) = Ro*diag(pH3)*pH2*pH1;
0380
0381 if SideInfo == 1
0382 if (i==1) & (j==n)
0383
0384 b3 = Q(II,J);
0385 EE = Q(I,J) + c3*A3h*Inv3*b3;
0386 for z=1:r
0387 for w=1:r
0388
0389
0390 F_rfc_z{z,w}(i,j) = 1/2*Ro(z)*EE(z,w);
0391
0392
0393 F_rfc_z{z,w}(j,i) = 1/2*Ro(z)*EE(z,w);
0394 end
0395 end
0396
0397 else
0398
0399 for z=1:r
0400 for w=1:r
0401
0402
0403 F_rfc_z{z,w}(i,j) = Roh(w)*pS3(w)*pS2(w,z)*pS1(z);
0404
0405
0406 F_rfc_z{z,w}(j,i) = Ro(w)*pH3(w)*pH2(w,z)*pH1(z);
0407
0408 end
0409 end
0410
0411 end
0412 elseif SideInfo ~= 0
0413 error(['SideInfo = ' num2str(SideInfo) ' not implemented']);
0414 end
0415
0416 end
0417 end
0418
0419
0420
0421
0422 F_rfc = c_m*F_rfc;
0423
0424 if SideInfo == 1
0425 for z=1:r
0426 for w=1:r
0427 F_rfc_z{z,w} = c_m*F_rfc_z{z,w};
0428 end
0429 end
0430 end
0431
0432 T(9,:)=clock;
0433
0434