0001
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 function [A,weights,dx,v,W,Nmatrix]=...
0037 estimate_dem_robust...
0038 (N,U,Np,Nr,Mc,poi,xa,weights,sigma_k,iter,type_robust,L1,...
0039 out_in,print_type,plot_type);
0040
0041
0042
0043
0044
0045 k=3;
0046
0047
0048
0049
0050 if print_type > 0
0051 tic
0052 start_time_sol = cputime;
0053 end
0054 A = spalloc(N,U,6*N);
0055 b = zeros(N,1);
0056 W = spalloc(N,N,N);
0057
0058 if print_type > 0
0059 truepos =sum(out_in .* (1-ceil(weights(1:Np)-0.5)));
0060 trueneg =sum((1-out_in) .* (ceil(weights(1:Np)-0.5)));
0061 falsepos=sum((1-out_in) .* (1-ceil(weights(1:Np)-0.5)));
0062 falseneg=sum(out_in .* (ceil(weights(1:Np)-0.5)));
0063
0064 Class_vor=[truepos,falseneg;falsepos,trueneg]
0065 end
0066
0067
0068 for p = 1:Np
0069 x = poi(p,1);
0070 y = poi(p,2);
0071 h = poi(p,3);
0072 w = sqrt(weights(p));
0073 s = poi(p,4);
0074 i = floor(x);
0075 j = floor(y);
0076 aa = x-i;
0077 bb = y-j;
0078
0079 A(p,i+j*Nr+1) = (1-aa)*(1-bb)/s*w;
0080 A(p,i+j*Nr+2) = aa *(1-bb)/s*w;
0081
0082 A(p,i+(j+1)*Nr+1) = (1-aa)* bb /s*w;
0083 A(p,i+(j+1)*Nr+2) = aa * bb /s*w;
0084 delta_l = h - ...
0085 ((1-aa)*(1-bb)*xa(i+j*Nr+1) +...
0086 aa *(1-bb)*xa(i+j*Nr+2) +...
0087 (1-aa)* bb *xa(i+(j+1)*Nr+1)+...
0088 aa * bb *xa(i+(j+1)*Nr+2));
0089 b(p) = delta_l /s*w;
0090 W(p,p)=1/s*w;
0091 end
0092
0093
0094
0095 if print_type > 0
0096 time_for_building_An=toc
0097 end
0098
0099
0100 if print_type > 0
0101 tic
0102 end
0103 Icc = Np;
0104 Irr = Np + (Nr-2)*Mc;
0105 Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);
0106
0107
0108
0109
0110 Jcc = [-1,0,+1];
0111 Vcc = [ 1 -2 1];
0112 wcc = 1/sigma_k;
0113
0114
0115 Jrr = [-Nr,0,+Nr];
0116 Vrr = [1 -2 1];
0117 wrr = 1/sigma_k;
0118
0119
0120
0121
0122 Jrc = [-Nr-1,-Nr,-1,0];
0123 Vrc = [1 -1 -1 1];
0124
0125 wrc = 1/sigma_k/sqrt(2);
0126
0127
0128 for j=1:Mc
0129 for i=2:Nr-1
0130 IU = i+(j-1)*Nr;
0131
0132 Icc = Icc+1;
0133 A(Icc,IU+Jcc)= Vcc * wcc * sqrt(weights(Icc));
0134 delta_cc = 0-xa(IU+Jcc)' * Vcc';
0135 b(Icc) = delta_cc * wcc * sqrt(weights(Icc));
0136 W(Icc,Icc) = wcc * sqrt(weights(Icc));
0137 end
0138 end
0139
0140 for j=2:Mc-1
0141 for i=1:Nr
0142 IU = i+(j-1)*Nr;
0143
0144 Irr = Irr+1;
0145 A(Irr,IU+Jrr)= Vrr * wrr * sqrt(weights(Irr));
0146 delta_rr = 0-xa(IU+Jrr)' * Vrr';
0147 b(Irr) = delta_rr * wrr * sqrt(weights(Irr));
0148 W(Irr,Irr)=wrr * sqrt(weights(Irr));
0149 end
0150 end
0151
0152 for j=2:Mc
0153 for i=2:Nr
0154 IU = i+(j-1)*Nr;
0155
0156 Irc = Irc+1;
0157 A(Irc,IU+Jrc) = Vrc * wrc * sqrt(weights(Irc));
0158 delta_rc = 0-xa(IU+Jrc)' * Vrc';
0159 b(Irc) = delta_rc * wrc * sqrt(weights(Irc));
0160 W(Irc,Irc)=wrc * sqrt(weights(Irc));
0161 end
0162 end
0163
0164
0165 if print_type > 0
0166 time_for_building_A=toc
0167 end
0168 if print_type > 1
0169 if U < 1000
0170 Ab = 2*[full(A)]
0171 null_A_prior = U-rank(full(A(U+1:end,:)));
0172 rankA = rank(full(A(U+1:end,:)));
0173 end
0174
0175
0176 end
0177
0178
0179
0180
0181 if print_type > 0
0182 tic
0183 end
0184 Nmatrix = A'*A;
0185 if iter==1 && plot_type > 1
0186 size_N_matrix = size(Nmatrix);
0187 non_zeros_N_matrix =nnz(Nmatrix);
0188 figure
0189 subplot(1,3,1)
0190 spy(Nmatrix)
0191 title('N')
0192 end
0193
0194 permx = symamd(Nmatrix);
0195 Aperm = A(:,permx);
0196
0197
0198 if iter==1 && plot_type > 1
0199 subplot(1,3,2)
0200 spy(Aperm'*Aperm)
0201 title('N, sorted')
0202 end
0203
0204
0205 [C,R] = qr(Aperm,b,0);
0206
0207 dxperm = R\C;
0208 Pm = speye(U);
0209 Pm = Pm(:,permx);
0210 dx = Pm*dxperm;
0211 if print_type > 0
0212 time_for_QR_sorted=toc
0213 end
0214
0215 if iter==1 && plot_type > 1
0216 subplot(1,3,3)
0217 spy(R)
0218 fill_in = nnz(R)-(nnz(Nmatrix)/2+U/2)
0219 percentage_nnz = nnz(R)/(U*(U+1)/2)
0220 relative_fill_in = fill_in/(nnz(Nmatrix)/2+U/2)
0221 title(strcat('reduced N sorted, fill in =',num2str(fill_in)))
0222 end
0223
0224
0225 if iter==1 && plot_type > 1
0226 figure
0227 subplot(1,2,1)
0228 spy(A)
0229 title('A')
0230 subplot(1,2,2)
0231 spy(Aperm)
0232 title('A sorted')
0233 end
0234
0235 dx_sorted = dx;
0236
0237 if U < 1600 && plot_type > 0
0238 Aqt = R\Aperm';
0239 figure
0240 spy(Aqt')
0241 title('R\A sorted')
0242 end
0243
0244 if print_type > 0
0245 time_for_solution=cputime-start_time_sol
0246 end
0247
0248
0249 v = A*dx-b;
0250
0251 eso = norm(v)/sqrt(N-U);
0252
0253 if type_robust > 0
0254 eso_p=median(abs(v(1:Np)))*1.48;
0255 vori = (v*sigma_k./sqrt(weights));
0256 vs=vori;
0257
0258 eso_k=median(abs(vs(Np+1:end)))*1.48;
0259 if print_type > 0
0260 est_s0_s0p_s0k=[eso,eso_p,eso_k]
0261 end
0262
0263 check_AtV_null=norm((A'*v));
0264 max_v_points = max(abs(v(1:Np)));
0265 max_v_curvat = max(abs(vs(Np+1,end)));
0266
0267 kp = k*eso_p;
0268 kp = k;
0269 kk = k*eso_k;
0270 if print_type > 0
0271 criticalvalue_kp_sigmap_kk_sigma_k=[kp,poi(1,4),kk,sigma_k]
0272 end
0273 end
0274
0275 if iter == 1 && type_robust > 0
0276 weights=ones(N,1);
0277 else
0278 if type_robust > 0 && L1
0279 if type_robust == 1 | type_robust == 3
0280 weights(Np+1:end) = min(1,1./(abs(vs(Np+1:end))/kk+eps));
0281 end
0282 if type_robust == 2 | type_robust ==3
0283 weights(1:Np) = min(1,1./(abs(v(1:Np))/kp+eps));
0284 end
0285 else
0286 if type_robust == 1 | type_robust == 3
0287 weights(Np+1:end) = exp(-abs(vs(Np+1:end)).^2/kk^2/2)+0.0001;
0288
0289 end
0290 if type_robust == 2 | type_robust ==3
0291 weights(1:Np) = exp(-abs(v(1:Np)).^2/kp^2/2)+0.0001;
0292 end
0293 end
0294 end
0295
0296 abs_dx = norm(dx);
0297 min_w_p=min(weights(1:Np));
0298 min_w_c=min(weights(Np+1:end));
0299
0300 if print_type > 1
0301 full(inv(W)*A)
0302 end
0303
0304
0305 if print_type > 0
0306 truepos =sum(out_in .* (1-ceil(weights(1:Np)-0.5)));
0307 trueneg =sum((1-out_in) .* (ceil(weights(1:Np)-0.5)));
0308 falsepos=sum((1-out_in) .* (1-ceil(weights(1:Np)-0.5)));
0309 falseneg=sum(out_in .* (ceil(weights(1:Np)-0.5)));
0310
0311 Class_nach=[truepos,falseneg;falsepos,trueneg]
0312 end
0313
0314 return