0001 function [x,z,TT] = smctpsim(P,F,T,init,whatOut)
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 TT(1,:) = clock;
0067
0068 ni = nargin;
0069 no = nargout;
0070 error(nargchk(3,5,ni));
0071
0072 Zstr = '123456789';
0073
0074 if ni < 4, init = []; end
0075 if ni < 5, whatOut = []; end
0076
0077 if isempty(init)
0078 init.x0 = [];
0079 init.z0 = [];
0080 end
0081 if ~isfield(init,'x0'), init.x0=[]; end
0082 if ~isfield(init,'z0'), init.z0=[]; end
0083
0084 if isempty(whatOut)
0085 whatOut = 'x';
0086 end
0087
0088 if isa(P,'double')
0089 Ptype = 'P';
0090 elseif isa(P,'struct')
0091 Ptype = 'struct';
0092 S = P;
0093 else
0094 error('P should be matrix or struct-array.');
0095 end
0096
0097 r = size(F,1);
0098
0099
0100 n = length(F{1,1});
0101
0102
0103
0104 if strcmp(Ptype,'struct')
0105 if isfield(S,'P')
0106 P=S.P;
0107 else
0108 P=[];
0109 end
0110 end
0111
0112 if ~isempty(P)
0113 sumP = sum(P');
0114 if sum(sumP == 1) ~= length(P)
0115 warning([': Rowsums of P not equal to 1. Renormalizing!']);
0116 for i = 1:length(P)
0117 P(i,:) = P(i,:)/sumP(i);
0118 end
0119 end
0120 end
0121
0122 TT(2,:) = clock;
0123
0124
0125
0126
0127 for i = 1:r
0128 QQ{i,1} = triu(F{i,1},1);
0129
0130 [I,J] = find(QQ{i,1}<0);
0131 if length(I) ~= 0
0132 warning(['Negative elements in Q' Zstr(i) '. Setting to zero!']);
0133 for k = 1:length(I)
0134 QQ{i,1}(I(k),J(k)) = 0;
0135 end
0136 end
0137
0138 sumQQi = sum(QQ{i,1}');
0139
0140 if sum(sumQQi == 1) ~= length(QQ{i,1})
0141
0142 for j = 1:n-1
0143 if sumQQi(j)~=0, QQ{i,1}(j,:) = QQ{i,1}(j,:)/sumQQi(j); end
0144 end
0145 end
0146 end
0147
0148 TT(3,:) = clock;
0149
0150
0151
0152
0153
0154
0155 for i = 1:r
0156 if isempty(F{i,2})
0157 QQ{i,2} = F{i,1}';
0158 else
0159 QQ{i,2} = F{i,2};
0160 end
0161
0162 QQ{i,2} = tril(QQ{i,2},-1);
0163
0164 [I,J] = find(QQ{i,2}<0);
0165 if length(I) ~= 0
0166 warning(['Negative elements in Qh' Zstr(i) '. Setting to zero!']);
0167 for k = 1:length(I)
0168 QQ{i,2}(I(k),J(k)) = 0;
0169 end
0170 end
0171
0172 sumQQi = sum(QQ{i,2}');
0173 if sum(sumQQi == 1) ~= length(QQ{i,2})
0174
0175 for j = 2:n
0176 if sumQQi(j)~=0, QQ{i,2}(j,:) = QQ{i,2}(j,:)/sumQQi(j); end
0177 end
0178 end
0179 end
0180
0181 TT(4,:) = clock;
0182
0183
0184
0185
0186
0187
0188
0189 [Q] = smctp2joint(P,QQ);
0190
0191
0192 [ro_min,ro_max] = mctp2stat(Q);
0193 ro = ro_min;
0194
0195
0196 e0 = rand(1,1);
0197 if isempty(init.z0) & isempty(init.x0)
0198 x0z0 = min(find( e0<=cumsum(ro) ));
0199 x0 = floor((x0z0+1)/r);
0200 z0 = mod(x0z0-1,r)+1;
0201 elseif isempty(init.x0)
0202 z0 = init.z0;
0203 rox0 = ro(z0:r:end);
0204 rox0 = rox0/sum(rox0);
0205 x0 = min(find( e0<=cumsum(rox0) ));
0206 elseif isempty(init.z0)
0207 x0 = init.x0;
0208 z0 = [];
0209 else
0210 x0 = init.x0;
0211 z0 = init.z0;
0212 end
0213
0214
0215 if strcmp(Ptype,'struct')
0216 if isfield(S,'P')
0217 S.P=P;
0218 end
0219 end
0220
0221
0222
0223
0224
0225
0226 switch whatOut
0227
0228 case {'x','x,RFM'}
0229 if strcmp(Ptype,'P')
0230 z = mcsim(P,T,z0);
0231 else
0232 z=zeros(T,1);
0233 z(1)=init_z(S,init);
0234 end
0235 e=rand(T,1);
0236 x=zeros(T,1);
0237
0238 case {'RFM'}
0239 if strcmp(Ptype,'P')
0240 z = mcsim(P,1,z0);
0241 else
0242 z=init_z(S,init);
0243 end
0244 e=rand(1,1);
0245
0246 end
0247
0248 x(1) = x0;
0249
0250 TT(5,:) = clock;
0251
0252
0253
0254
0255
0256
0257
0258 cumsumP = cumsum(P,2);
0259 ones_n = ones(n,1);
0260 for i = 1:r
0261 cumsumQQ{i,1} = cumsum(QQ{i,1},2);
0262 cumsumQQ{i,2} = cumsum(QQ{i,2},2);
0263 cumsumQQ{i,1}(:,end) = ones_n;
0264 cumsumQQ{i,2}(:,end) = ones_n;
0265 end
0266
0267
0268
0269 switch whatOut
0270
0271 case {'x','x,RFM'}
0272 switch Ptype
0273 case 'P'
0274 for k=2:T
0275 if rem(k,2) == 0
0276 x(k) = sum( cumsumQQ{z(k),1}(x(k-1),:) <= e(k) ) + 1;
0277 else
0278 x(k) = sum( cumsumQQ{z(k),2}(x(k-1),:) <= e(k) ) + 1;
0279 end
0280 end
0281
0282 case 'struct'
0283 t0 = simT(S,z(1));
0284 for k=2:T
0285 if t0 > 1
0286 z(k) = z(k-1);
0287 t0=t0-1;
0288 else
0289 z(k-1:k) = mcsim(S.P,2,z(k-1));
0290 t0 = simT(S,z(k));
0291 end
0292
0293 if rem(k,2) == 0
0294 x(k) = sum( cumsumQQ{z(k),1}(x(k-1),:) <= e(k) ) + 1;
0295 else
0296 x(k) = sum( cumsumQQ{z(k),2}(x(k-1),:) <= e(k) ) + 1;
0297 end
0298 end
0299 end
0300
0301
0302 case 'test'
0303 [RFM0,res,nres] = dtp2rfm_init(n);
0304 [RFM0,res,nres] = dtp2rfm1(x(1),RFM0,res,nres);
0305 for k=2:T
0306 e=rand(2,1);
0307
0308
0309 z(k) = sum( cumsumP(z(k-1),:) <= e(1) ) + 1;
0310
0311
0312 if rem(k,2) == 0
0313 x(k) = sum( cumsumQQ{z(k),1}(x(k-1),:) <= e(2) ) + 1;
0314 else
0315 x(k) = sum( cumsumQQ{z(k),2}(x(k-1),:) <= e(2) ) + 1;
0316 end
0317 [RFM0,res,nres] = dtp2rfm1(x(k),RFM0,res,nres);
0318 end
0319 [RFM] = dtp2rfm_collect(RFM0,res,nres);
0320
0321 case 'RFM'
0322 [RFM0,res,nres] = dtp2rfm_init(n);
0323 [RFM0,res,nres] = dtp2rfm1(x,RFM0,res,nres);
0324 for k=2:T
0325 e=rand(2,1);
0326
0327
0328 z = sum( cumsumP(z,:) <= e(1) ) + 1;
0329
0330
0331 if rem(k,2) == 0
0332 x = sum( cumsumQQ{z,1}(x,:) <= e(2) ) + 1;
0333 else
0334 x = sum( cumsumQQ{z,2}(x,:) <= e(2) ) + 1;
0335 end
0336 [RFM0,res,nres] = dtp2rfm1(x,RFM0,res,nres);
0337 end
0338 [RFM] = dtp2rfm_collect(RFM0,res,nres);
0339
0340 end
0341
0342 TT(6,:) = clock;
0343
0344
0345
0346 switch whatOut
0347 case 'x,RFM'
0348 RFM = dtp2rfm(x,n);
0349 TT = RFM;
0350 case 'RFM'
0351 x = RFM;
0352 end
0353
0354
0355 function z=init_z(S,init)
0356
0357
0358
0359 if isempty(init.z0)
0360 r=length(S.P);
0361 m = zeros(1,r);
0362 switch S.distr
0363 case {'exp','const'},
0364 for i=1:r, m(i)=S.Tpar{i}; end
0365 case 'phasetype',
0366 for i=1:r, m(i)=sum(S.Tpar{i}); end
0367 end
0368 statP=mc2stat(S.P);
0369 ro = statP.*m;
0370 z = sum( cumsum(ro) <= rand ) + 1;
0371 else
0372 z=init.z0;
0373 end
0374
0375
0376 function t=simT(S,z)
0377
0378
0379
0380 if iscell(S.distr)
0381 distr = S.distr{z};
0382 else
0383 distr = S.distr;
0384 end
0385
0386 switch distr
0387
0388 case 'exp'
0389 t = simFS(S.Tpar{z});
0390 case 'phasetype'
0391 t = simPhaseType(S.Tpar{z});
0392 case 'const'
0393 t = S.Tpar{z};
0394 case 'other'
0395 error('other: Not yet implemented.');
0396 end
0397
0398
0399
0400 function t=simFS(m)
0401
0402
0403 u = rand(size(m));
0404 p = 1./m;
0405 t = floor(log(u) ./ log(1 - p)) + 1;
0406
0407
0408
0409 function t=simExp(m)
0410
0411 t = - m .* log(rand(size(m)));
0412
0413
0414
0415 function t=simPhaseType(m)
0416
0417 t = sum(simFS(m));
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429 function [RFM0,res,nres] = dtp2rfm_init(n)
0430
0431
0432
0433 RFM0 = zeros(n);
0434 res = zeros(2*n+1,1);
0435 nres = 0;
0436
0437
0438 function [RFM0,res,nres] = dtp2rfm1(x,RFM0,res,nres)
0439
0440
0441
0442
0443
0444
0445
0446 nres = nres+1;
0447 res(nres) = x;
0448 cycleFound = 1;
0449 while cycleFound==1 & nres>=4
0450 A = sort([res(nres-1) res(nres-2)]);
0451 B = sort([res(nres) res(nres-3)]);
0452 if A(1) >= B(1) & A(2) <= B(2)
0453 RFM0(res(nres-2),res(nres-1)) = RFM0(res(nres-2),res(nres-1)) + 1;
0454 res(nres-2) = res(nres);
0455 nres = nres-2;
0456 else
0457 cycleFound = 0;
0458 end
0459 end
0460
0461
0462
0463
0464 function [RFM] = dtp2rfm_collect(RFM0,res,nres)
0465
0466
0467
0468
0469
0470
0471 RFM = RFM0;
0472 for i=1:2:nres-1
0473 RFM(res(i),res(i+1)) = RFM(res(i),res(i+1)) + 1;
0474 end
0475
0476
0477
0478 RFM = RFM+RFM';
0479 RFM = triu(RFM);
0480
0481