0001 function [Fest,Pout,Fextreme,Fsmooth,Fest0] = rfmextrapolate2(F,Pin,plotflag)
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 ni = nargin;
0068 no = nargout;
0069 error(nargchk(1,3,ni));
0070
0071 if ni<2, Pin=[]; end
0072 if ni<3, plotflag=[]; end
0073
0074 if isempty(plotflag)
0075 plotflag=0;
0076 end
0077
0078 n = length(F);
0079
0080
0081
0082
0083
0084
0085 Pout.param = [1 n n];
0086
0087 Pout.method = 'gpd,ml';
0088 Pout.u_lev = [];
0089 Pout.LCfrac = [];
0090
0091 Pout.h = -1;
0092
0093 Pout.Lim = [];
0094 Pout.LimRelDam = [];
0095 Pout.beta = 7;
0096
0097
0098 Pout.PL = [10:20:90 99 99.9 99.99 99.999];
0099
0100
0101 if ~isempty(Pin)
0102 Fname = fieldnames(Pin);
0103 for i = 1:length(Fname)
0104 Pout = setfield(Pout,Fname{i},getfield(Pin,Fname{i}));
0105 end
0106 end
0107
0108
0109
0110
0111
0112 if isempty(Pout.u_lev) & isempty(Pout.LCfrac)
0113 Pout.LCfrac = 0.05;
0114 end
0115 if isempty(Pout.Lim) & isempty(Pout.LimRelDam)
0116 Pout.LimRelDam = 0.95;
0117 end
0118
0119
0120
0121 param = Pout.param;
0122 paramD = [1 n n];
0123 beta = Pout.beta;
0124 PL = Pout.PL;
0125
0126 uD = levels(paramD);
0127 u = levels(param);
0128 du=u(2)-u(1);
0129
0130
0131
0132
0133
0134
0135 LCfrac = Pout.LCfrac;
0136 u_lev = Pout.u_lev;
0137
0138 if plotflag >=2, figure, end
0139
0140 if isempty(u_lev)
0141 slut=0;
0142 F1=F; i_lowPrev=0; i_highPrev=0;
0143 while ~slut
0144 lc = cmat2lc(paramD,F1);
0145
0146 methodFindLev=3;
0147 if methodFindLev == 1
0148 [M,Mi] = max(lc(:,2)); Mi = Mi(1);
0149 dlc = [diff([0; lc(1:Mi-1,2)]); 0; flipud(diff([0; flipud(lc(Mi+1:end,2))]))];
0150
0151 Nlc = sum(dlc~=0);
0152 Nextr = floor(0.2*Nlc);
0153 if Nextr<3, Nextr=3; end
0154
0155 I = find(dlc~=0);
0156 i_low0 = I(Nextr+1);
0157 i_high0 = I(end-Nextr);
0158
0159 elseif methodFindLev == 2
0160
0161 [M,Mi] = max(lc(:,2)); Mi = Mi(1);
0162 dlc1 = diff([0; lc(1:Mi-1,2)]);
0163 dlc2 = diff([0; flipud(lc(Mi+1:end,2))]);
0164
0165 Nlc1 = sum(dlc1~=0);
0166 Nlc2 = sum(dlc2~=0);
0167 Nextr1 = max(ceil(0.4*Nlc1),3);
0168 Nextr2 = max(ceil(0.4*Nlc2),3);
0169
0170
0171 I1 = find(dlc1~=0);
0172 I2 = find(dlc2~=0);
0173 i_low0 = I1(Nextr1+1);
0174 i_high0 = n-I2(Nextr2+1)+1;
0175
0176 I = [I1; n-flipud(I2)+1];
0177
0178 elseif methodFindLev == 3
0179
0180 [M,Mi] = max(lc(:,2)); M = M(1); Mi=Mi(1);
0181
0182 I = find(lc(:,2)>LCfrac*M);
0183 i_low0 = I(1);
0184 i_high0 = I(end);
0185
0186 end
0187
0188 [M,Mi] = max(lc(:,2)); M = M(1); Mi=Mi(1);
0189 dlc = [diff([0; lc(1:Mi-1,2)]); 0; flipud(diff([0; flipud(lc(Mi+1:end,2))]))];
0190 Nlc = sum(dlc~=0);
0191 Nextr = 3;
0192 I = find(dlc~=0);
0193 i_low = max(i_low0,I(Nextr+1));
0194 i_high = min(i_high0,I(end-Nextr));
0195
0196
0197 F1 = F;
0198 for i = i_high+1:n
0199 F1(i,:) = 0;
0200 end
0201 for j = 1:i_low-1
0202 F1(:,j) = 0;
0203 end
0204
0205 slut = (i_high==i_highPrev) & (i_low==i_lowPrev);
0206 i_highPrev=i_high; i_lowPrev=i_low;
0207
0208 end
0209 else
0210 i_low = u_lev(1); i_high = u_lev(2);
0211
0212 F1 = F;
0213 for i = i_high+1:n
0214 F1(i,:) = 0;
0215 end
0216 for j = 1:i_low-1
0217 F1(:,j) = 0;
0218 end
0219 lc = cmat2lc(paramD,F);
0220 end
0221
0222 if plotflag >=2
0223 [M,Mi] = max(lc(:,2)); M = M(1); Mi=Mi(1);
0224 lc1 =lc; lc1(lc(:,2)==0,2) = 0.1;
0225 stairs(u(lc1(1:Mi,1)),lc1(1:Mi,2)), hold on
0226 stairs(u(lc1(Mi+1:end,1)-1),lc1(Mi+1:end,2)), hold on
0227
0228
0229
0230
0231 v = axis; axis([u([1 n]) 0.9 v(4)]);
0232
0233 plot(u([i_low i_high]),0.9*ones(1,2),'r.'),
0234 plot(u([i_low i_low]),[0.9 v(4)],'r'),
0235 plot(u([i_high i_high]),[0.9 v(4)],'r'),
0236 if ~isempty(LCfrac), plot(u([1 n]),LCfrac*[M M],'r'), end
0237 hold off
0238 set(gca,'Yscale','log')
0239 title('Choice of threshold for extrapolation')
0240 xlabel('level'), ylabel('Number of crossings')
0241 end
0242
0243
0244
0245
0246 u_lev = [i_low i_high];
0247 if plotflag>=2
0248 plotflag2 = 1; figure
0249 else
0250 plotflag2 = 0;
0251 end
0252
0253 method = Pout.method;
0254 [lcEst,Est] = cmat2extralc(paramD,F,u_lev,method,plotflag2);
0255
0256 if plotflag >= 2
0257 figure
0258 lcH =lcEst.High; lcH(lcH(:,2)==0,2) = 1e-100;
0259 lcL =lcEst.Low; lcL(lcL(:,2)==0,2) = 1e-100;
0260 plot(u(lc(:,1)),lc(:,2),'-'), hold on
0261 plot(u(lcH(:,1)),lcH(:,2),'r')
0262 plot(u(lcL(:,1)),lcL(:,2),'r')
0263 hold off, grid on
0264 set(gca,'YScale','log')
0265 Mlc = 10^(ceil(log10(max(lc(:,2)))));
0266 v=axis; axis([u([1 n]) Mlc/1e8 Mlc])
0267 title('Extrapolated level crossings')
0268 xlabel('level, u')
0269 ylabel('Level crossing intensity, \mu(u)')
0270 end
0271 Pout.u_lev = u_lev;
0272
0273
0274
0275
0276
0277 lc1 = [0 0; lcEst.lc; n+1 0];
0278
0279 Ga1 = lc2rfmextreme(lc1);
0280 Ga=Ga1(2:n+1,2:n+1);
0281 Ga(isnan(Ga)) = 0;
0282
0283 if plotflag >= 2
0284 figure
0285 [qlGa PL] = qlevels(F,PL,u,u);
0286
0287 contour(u,fliplr(u),flipud(F'),qlGa,'b'), hold on
0288 contour(u,fliplr(u),flipud(Ga'),qlGa,'r'), hold off
0289 title('Extreme RFM (red), compared with observed RFM (blue)')
0290 xlabel('min')
0291 ylabel('Max')
0292 v = axis; axis([min(u) max(u) min(u) max(u)]);
0293 axis('square')
0294 end
0295
0296
0297
0298
0299
0300 h=Pout.h;
0301
0302 if h~=0
0303
0304 i=0;
0305 while sum(diag(F,i)) == 0
0306 i=i+1;
0307 end
0308 NOsubzero = i-1;
0309
0310 if h==-1
0311 h_norm = smoothcmat_hnorm(F,NOsubzero);
0312
0313
0314
0315 h=0.5*h_norm;
0316 end
0317
0318 [Gs] = smoothcmat(F,1,h,NOsubzero);
0319
0320 else
0321 Gs=F;
0322 end
0323
0324 Pout.h = h;
0325
0326 if plotflag >= 2
0327 figure
0328 [qlGs PL] = qlevels(Gs,PL,u,u);
0329
0330 contour(u,fliplr(u),flipud(Ga'),qlGs,'r'), hold on
0331 contour(u,fliplr(u),flipud(Gs'),qlGs,'b'), hold off
0332 title('Smooth RFM (blue), compared with Extreme RFM (red)')
0333 xlabel('min')
0334 ylabel('Max')
0335 v = axis; axis([min(u) max(u) min(u) max(u)]);
0336 axis('square')
0337 end
0338
0339
0340
0341
0342
0343 Dnorm = cmat2dam(paramD,F,beta)*100;
0344 ampF = cmat2amp(paramD,cmat2dmat(paramD,F/Dnorm,beta));
0345 ampGa = cmat2amp(paramD,cmat2dmat(paramD,Ga/Dnorm,beta));
0346
0347 ptyp=2;
0348 Drel = zeros(n,1);
0349 IampF = (ampF(:,2)==0);
0350 IampGa = (ampGa(:,2)==0);
0351 I = ~IampF & ~IampGa;
0352 if ptyp < 3
0353 Drel(I)=ampGa(I,2)./ampF(I,2);
0354 Drel(IampF)=1e10;
0355 else
0356 Drel(I)=ampF(I,2)./ampGa(I,2);
0357 Drel(IampGa)=1e10;
0358 end
0359
0360 Drel(IampF & IampGa)=0;
0361
0362
0363 if plotflag >= 2
0364 figure
0365 subplot(2,1,1),plot(ampF(:,1)*du,ampF(:,2)),hold on
0366 plot(ampGa(:,1)*du,ampGa(:,2),'r'),hold off
0367 v = axis; axis([0 n/2*du v(3:4)])
0368 title('Damage per amplitude'), xlabel('amplitude'), ylabel('D0, Dextreme')
0369 subplot(2,1,2),plot(ampF(:,1)*du,Drel), hold on
0370 if ptyp<3
0371 plot([0 n/2*du],0.95*[1 1],'r--'), hold off
0372
0373 v = axis; axis([0 n/2*du 0 5])
0374 title('Relative damage per amplitude'), xlabel('amplitude'), ylabel('Drel = Dextreme / D0')
0375 if ptyp==2, set(gca,'YDir','reverse'), end
0376 elseif ptyp==3
0377 plot([0 n/2*du],1/0.95*[1 1],'r--'), hold off
0378 v = axis; axis([0 n/2*du 0 5])
0379 title('Relative damage per amplitude'), xlabel('amplitude'), ylabel('Drel = D0 / Dextreme')
0380 end
0381
0382 end
0383
0384 Lim = Pout.Lim;
0385
0386 if ~isfield(Lim,'range')
0387
0388
0389
0390 LimRelDam = Pout.LimRelDam;
0391 I = find((ampF(:,2)~=0) & (ampGa(:,2)~=0));
0392 k = min(find(Drel(I)>LimRelDam));
0393 if ~isempty(k)
0394 Lim.range = I(k);
0395 else
0396 Lim.range=n;
0397 end
0398 end
0399
0400 [dummy,imax] = max(lcEst.lc(:,2));
0401 maxLC = imax(1);
0402 gap=floor(0.02*n)+1;
0403 if ~isfield(Lim,'min')
0404 Lim.min = maxLC-gap;
0405 end
0406 if ~isfield(Lim,'max')
0407 Lim.max = maxLC+gap;
0408 end
0409
0410 [GaC0,LimOut] = cmatcombine(Ga,Gs,Lim);
0411
0412 Pout.Lim = Lim;
0413
0414
0415
0416
0417 f_max = sum(Ga);
0418 jj=min(find([0 fliplr(f_max)]>0));
0419 jmax = n-jj+2;
0420
0421 f_min = sum(Ga');
0422 ii=min(find([0 f_min]>0));
0423 imin = ii-1;
0424
0425 J=meshgrid(1:n);
0426 I=J';
0427 K1 = (I < imin) | (J > jmax);
0428 GaC = GaC0;
0429 GaC(K1) = 0;
0430
0431 Pout.imin = imin;
0432 Pout.jmax = jmax;
0433
0434 if plotflag >= 1
0435 if plotflag>=2, figure, end
0436 [qlGaC PL] = qlevels(GaC,PL,u,u);
0437
0438 contour(u,fliplr(u),flipud(F'),qlGaC,'b'), hold on
0439 contour(u,fliplr(u),flipud(GaC'),qlGaC,'r'),
0440
0441 plot(u([1 n-LimOut.range+1]),u([LimOut.range n]),'g')
0442 plot(u([1 n]),u([LimOut.max LimOut.max]),'g')
0443 plot(u([LimOut.min LimOut.min]),u([1 n]),'g'),hold off
0444 title('Limiting RFM (red), compared with observed RFM (blue)')
0445 xlabel('min')
0446 ylabel('Max')
0447 v = axis; axis([min(u) max(u) min(u) max(u)]);
0448 axis('square')
0449 end
0450
0451
0452
0453 if plotflag >= 2
0454 figure
0455
0456 DmatGs = cmat2dmat(paramD,Gs,beta);
0457 DmatGa = cmat2dmat(paramD,Ga,beta);
0458 DmatGaC = cmat2dmat(paramD,GaC,beta);
0459 DmatF = cmat2dmat(paramD,F,beta);
0460
0461 DamF = sum(sum(DmatF));
0462 DrelGa = sum(sum(DmatGa))/DamF;
0463 DrelGs = sum(sum(DmatGs))/DamF;
0464 DrelGaC = sum(sum(DmatGaC))/DamF;
0465
0466 subplot(2,2,1), cmatplot(u,u,DmatF,13), title(['RFM, Drel = 1'])
0467 subplot(2,2,2), cmatplot(u,u,DmatGa,13), title(['Extreme RFM, Drel = ' num2str(DrelGa)])
0468 subplot(2,2,3), cmatplot(u,u,DmatGs,13), title(['Smoothed RFM, Drel = ' num2str(DrelGs)])
0469 subplot(2,2,4), cmatplot(u,u,DmatGaC,13), title(['Extrapolated RFM, Drel = ' num2str(DrelGaC)])
0470
0471 end
0472
0473
0474
0475
0476
0477 Fest = GaC;
0478
0479 if no>2
0480 Fest0 = GaC0;
0481 Fextreme = Ga;
0482 Fsmooth = Gs;
0483 end
0484
0485