cycling/squish4.m

   1 function tign_new = squish4(ps,gq,da)
   2 % ps is struct with paths, red, graph,distances, etc.
   3 %  gq  - use ground detetcions if gq == 1
   4 %  da - use in data assimilation mode if da == 1
   5 % ps = graph_dets(w)
   6
   7 if ~exist('a','var')
   8     a = 300;
   9     b = 1/3;
  10 end
  11
  12 %set flags for data assimilation mode
  13 if da == 1
  14     % controls how ground detections are handled
  15     use_beta_likes = 1;
  16     % control how active fire detections handled
  17     use_alpha_likes = 1;
  18 else
  19     use_beta_likes = 0;
  20     use_alpha_likes = 0;
  21 end
  22
  23 %combine l2 time with detection points fixed to grid
  24 pts = ps.grid_pts;
  25 pts(:,3) = ps.points(:,3);
  26 %adjust time of ignition
  27 %pts(1,3) = pts(1,3)-1/2;
  28 %forecast
  29 tign = ps.red.tign;
  30 %experiments
  31 % sat_max = max(ps.points(:,3));
  32 % sat_msk = tign > sat_max;
  33 % tign(sat_msk) = sat_max;
  34
  35 % t_max = max(tign(:))-.1;
  36 % tign(tign>=t_max) =t_max+1;
  37 t0 = min(tign(:));
  38 [m,n]=size(ps.red.tign);
  39
  40 % flat_tign = input_num('Use flat start? 1 = yes',0)
  41 % if flat_tign
  42 %     tign = ps.red.end_datenum*ones(m,n);
  43 % end
  44
  45 %%% data likelihood spline %%%
  46 %set up data likelihood spline
  47 [p_like_spline,~,~] = make_spline(100,1000);
  48 %test
  49 % t = 4*(-10:20)*3600;
  50 % ts = p_like_spline(t);
  51 % figure,plot(t,exp(ts))
  52 % hold on,plot(t,1-exp(ts))
  53
  54 %ground detection usage
  55 if ~exist('gq','var')
  56     gq = input_num('use ground detections? yes = [1]',0,1);
  57 end
  58 if gq
  59     smooth_ground = 1;%max(m,n)/50;
  60     fprintf('Collecting ground detection data\n')
  61
  62 %     k = boundary(ps.points(:,2),ps.points(:,1),1);
  63 %     lon_perim = ps.points(k,2);
  64 %     lat_perim = ps.points(k,1);
  65 %     in = inpolygon(ps.red.fxlat,ps.red.fxlong,lat_perim,lon_perim);
  66 %     infire = in;
  67
  68     %grounds = ground_detects(ps.red);
  69     %[jgrid,igrid]=meshgrid([1:length(ps.red.jspan)]',[1:length(ps.red.ispan)]');
  70     %[jgrid,igrid]=meshgrid([1:n]',[1:m]');
  71     k = boundary(ps.points(:,2),ps.points(:,1),1);
  72     lon_perim = ps.points(k,2);
  73     lat_perim = ps.points(k,1);
  74     in = inpolygon(ps.red.fxlong,ps.red.fxlat,lon_perim,lat_perim);
  75     infire = in;
  76     %make polygon around detections
  77     %infire = inpolygon(grounds.land(:,4),grounds.land(:,3),pts(:,2),pts(:,1));
  78     %infire = inpolygon(grounds.land(:,4),grounds.land(:,3),grounds.land(infire,4),grounds.land(infire,3));
  79     %infire = inpolygon(igrid(:),jgrid(:),ps.idx(:,1),ps.idx(:,2));
  80     %remove holes in mask
  81     %infire = inpolygon(igrid(:),jgrid(:),igrid(infire),jgrid(infire));
  82 end
  83 %try
  84
  85 tmax=max(tign(:));
  86 tign_ground = tign;
  87 %tign_flat is a tign with time equal to simulation end
  88 tign_flat = ps.red.end_datenum*ones(size(ps.red.tign));
  89 %figure(159),hold on,scatter(pts(:,2),pts(:,1),'*r')
  90 time_err = 0.0;
  91 g_diff = (tign_flat-tign+time_err)*24*3600;
  92 %
  93 g_likes = p_like_spline(g_diff);
  94 beta = 1/100;
  95 beta_vect = 1-exp(g_likes);
  96 %%% ground detection likelihood
  97 % g_times = zeros(pts_length,1);
  98 % for i = 1:pts_length
  99 %     g_times(i) = ps.red.tign(ps.idx(i,1),ps.idx(i,2));
 100 % end
 101 ground_steps = 3;
 102 if gq
 103     data_area = sum(infire(:));
 104     for i = 1:ground_steps
 105         fire_mask = tign_ground < tmax-0.1;
 106         fire_area = sum(fire_mask(:));
 107         out_side = fire_mask(:)-infire(:);
 108         out_side(out_side<0)=0;
 109         out_sum(i) = sum(out_side(:));
 110         out_fraction = out_sum(i)/data_area;
 111         if i > 2 && out_sum(i) == out_sum(i-1)
 112             smooth_ground = 0.9*smooth_ground;
 113             fprintf('decreasing ground smoothing \n')
 114         end
 115         if out_fraction < 0.1
 116             break
 117         end
 118         fprintf('Fire area: %f pixel_out/data area: %f  pixels_out %f\n',fire_area,out_sum(i)/data_area,out_sum(i));
 119         if use_beta_likes == 1
 120             tign_ground(~infire) = beta_vect(~infire).*tign_ground(~infire)+(1-beta_vect(~infire)).*tign_flat(~infire);
 121         else
 122             tign_ground(~infire) = beta*tign_ground(~infire)+(1-beta)*tign_flat(~infire);
 123         end
 124         %flatten the top
 125         t_mask = tign_ground > tmax;
 126         tign_ground(t_mask) = tmax;
 127         %tign_temp = imgaussfilt(tign_ground,smooth_ground );
 128         tign_temp = smooth_up(tign_ground,a,b);
 129         tign_ground(~infire) = tign_temp(~infire);
 130 %         a = 1/10;%1/2-1/(2*i);
 131 %         tign_ground(infire) = a*tign_ground(infire)+(1-a)*tign_temp(infire);
 132 %         figure(73),scatter3(ps.red.fxlong(~infire),ps.red.fxlat(~infire),tign_ground(~infire))
 133 %         pause(1/2)
 134 %         figure(159)
 135 %         contourf(ps.red.fxlong,ps.red.fxlat,tign_ground,20,'k'),hold on
 136 %         scatter(pts(1:10:end,2),pts(1:10:end,1),'*r'),hold off
 137 %         t_str = sprintf('Perimeter Shrinking \n Iteration %d',i);
 138 %         save_str = sprintf('perim_shrink_%d',i);
 139 %        xlabel('Lon [degree]'),ylabel('Lat [degre]'),title(t_str)
 140 %         savefig(save_str);
 141 %         saveas(gcf,[save_str '.png']);
 142         %figure(160),mesh(ps.red.fxlong,ps.red.fxlat,tign_ground);
 143         %pause(.5)
 144         %t_min(i) = min(tign_ground(:));
 145 %         close 159
 146     end
 147 end
 148 %plot(t_min)
 149 % figure,scatter(igrid(infire),jgrid(infire))
 150 % time_mask = grounds.land(:,5) > max(ps.points(:,3));
 151 % ground_mask = logical(time_mask.*~infire);
 152 tign_new = tign_ground;
 153 %%%plot and save original perimeter
 154 % figure(159),hold on
 155 % scatter(pts(:,2),pts(:,1),'*r')
 156 % contourf(ps.red.fxlong,ps.red.fxlat,ps.red.tign,20,'k'),hold off
 157 % t_str = sprintf('Perimeter Shrinking \n Iteration %d',0);
 158 % save_str = sprintf('perim_shrink_%d',0);
 159 % xlabel('Lon [degree]'),ylabel('Lat [degre]'),title(t_str)
 160 % xlim(xl);
 161 % ylim(yl);
 162 % savefig(save_str);
 163 % saveas(gcf,[save_str '.png']);
 164
 165 %%%%%%% end ground detection work %%%%%%
 166
 167
 168 idx = ps.idx;
 169 fig_num = 179;
 170 pts_length = length(ps.grid_pts);
 171 %max time to look at detection data
 172 max_l2_time = max(max(pts(:,3)));
 173
 174
 175
 176
 177 %tign_flat is constant fire arrival time
 178 %tign_flat=ps.red.end_datenum*ones(size(tign));
 179 title_str = 'Analysis';
 180
 181 %make matrix of forecast tign time for the points in the graph
 182 t_times = zeros(pts_length,1);
 183 for i = 1:pts_length
 184     t_times(i) = tign(idx(i,1),idx(i,2));
 185 end
 186
 187
 188
 189 %data likelikehoods
 190 time_diff = (pts(:,3)-t_times)*24*3600;
 191 likes = p_like_spline(time_diff);
 192 alpha_vect = exp(likes);
 193
 194 %mask = ones(size(tign));
 195 %change back to z = tign
 196 % if exist('ps.tn','var')
 197 %     figure(68),contour3(ps.red.fxlong,ps.red.fxlat,ps.tn-ps.red.start_datenum,20)
 198 % else
 199 %     figure(68),contour3(ps.red.fxlong,ps.red.fxlat,tign-ps.red.start_datenum,20)
 200 % end
 201 % hold on
 202 %lons = ps.red.fxlong;
 203 %lats = ps.red.fxlat;
 204
 205 %just squish the ignition point
 206 %bigger k ==> more smoothing
 207
 208 norms=[];
 209
 210 %random multiplier, increase for larger grids
 211 %perturbs points on path in x-y plane
 212 rm = 0;
 213 % random multiplier, keep the same
 214 % perturbs points downward in time to
 215 rt = 0;
 216 % weight for tign_new
 217
 218 %alhpa blends  estimate of tign at a point with old estimate
 219 % new_estimate = alph*current_setimate + (1-alpha)*old_estimate
 220 % for data assimilation, alpha will be computed from exp(likelihood)
 221 % alpha = 0.5;
 222 %constant for smooth in rlx_shp
 223 %alpha_2 = 0.7; %smaller alph_2 ==> smoother
 224 %number of loops to run
 225 smoothings = 20;
 226 if da ==1
 227     smoothings = 20;
 228 end
 229 for k = 1:smoothings
 230 %     figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new)
 231 %     title(title_str)
 232     %hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3),'*r'),hold off
 233     %pause(3/k)
 234     %plotting a slice
 235     %figure(fig_num+17),plot(tign_new(:,202));
 236     for i = 1:length(ps.paths)
 237         p = ps.paths(i).p;
 238         %figure(73),hold on
 239         %plot3(pts(p,2),pts(p,1),pts(p,3)-ps.red.start_datenum,'r')
 240         %hold off
 241         %plot3(pts(p,2),pts(p,1),tign(idx(p,1),idx(p,2))-ps.red.start_datenum,'g')
 242         for j = 1:length(p) %p(j) is the current detection vertex
 243             tign_old = tign_new;
 244             %mesh indices for path points, perturbed
 245             p_i = idx(p(j),1);%+rm*round(randn);
 246             p_j = idx(p(j),2);%+rm*round(randn);
 247             %%% make mean of old and new, in small block around path point
 248             %alpha is now the data likelikehood
 249             if use_alpha_likes == 1
 250                 alpha = alpha_vect(p(j));
 251             else
 252                 alpha = 0.0;
 253             end
 254             % adjust tign as weighted average of tign_new and time of detection
 255             tign_new(p_i-round(rm*rand):p_i+round(rm*rand),p_j-round(rm*rand):p_j+round(rm*rand)) = alpha*tign_new(p_i,p_j) + (1-alpha)*pts(p(j),3)-rt;
 256             %
 257             %tign_flat(p_i-round(rm*rand):p_i+round(rm*rand),p_j-round(rm*rand):p_j+round(rm*rand)) = 0.5*(tign_new(p_i,p_j) + pts(p(j),3)-rt*rand);
 258
 259             %no need to interpolate - already done
 260             % interpolate new point between adjacent path detections
 261 %             if j > 1
 262 %                 %weighted average will move new point close to that with
 263 %                 %higher FRP
 264 %                 %                 frp1 = ps.points(p(j-1),5);
 265 %                 %                 frp2 = ps.points(p(j),5);
 266 %                 %                 w1 = frp1/(frp1+frp2);
 267 %                 %                 w2 = frp2/(frp1+frp2);
 268 %                 w1 = 0.5;
 269 %                 w2 = 0.5;
 270 %                 %assign tign for all in small, block around midpoint
 271 %                 %             frp1 = ps.points(p(j-1),5);
 272 %                 %             frp2 = ps.points(p(j),5);
 273 %                 %             w1 = frp1/(frp1+frp2);
 274 %                 %             w2 = frp2/(frp1+frp2);
 275 %                 new_lat = w1*pts(p(j-1),1)+w2*pts(p(j),1);
 276 %                 new_lon = w1*pts(p(j-1),2)+w2*pts(p(j),2);
 277 %                 new_t = w1*pts(p(j-1),3)+w2*pts(p(j),3);
 278 %                 [new_i,new_j,new_lat,new_lon]= fixpt(ps.red,[new_lat,new_lon]);
 279 %                 tign_new(new_i-round(rm*rand):new_i+round(rm*rand),new_j-round(rm*rand):new_j+round(rm*rand)) = new_t-rt*rand;
 280 %             end %interpolate a point block
 281         end
 282     end
 283     %patch sets size of a local averaging in experimental process
 284     %patch = max(1,round(sqrt(smoothings-k)));
 285     %patch = 4*ceil(max(m,n)/100);
 286     %fprintf('Using patch size %d in rlx_shp function \n',patch);
 287     %patch = 2;
 288
 289     %%%%% make use of ground detections %%%%
 290     %tign_new(~infire) = beta*tign_new(~infire)+(1-beta)*tign_flat(~infire);
 291
 292     %smooth the tign
 293     %tign_new(tign_new < t0) = t0;
 294     tign_new = smooth_up(tign_new,a,b);
 295 %     if k < smoothings/5
 296 %         tign_new = imgaussfilt(tign_new,1);
 297 %         tign_new = rlx_shp(tign_new,alpha_2,patch);
 298 %     end
 299 %   tign_flat = rlx_shp(tign_flat,alpha_2,patch);
 300
 301     %collect information about tign at the detection points
 302     for i = 1:pts_length
 303         t_times(i) = tign_new(ps.idx(i,1),ps.idx(i,2));
 304         %flat_times(i) = tign_flat(ps.idx(i,1),ps.idx(i,2));
 305     end
 306     norms(k,1) = norm(t_times-ps.points(:,3),2);
 307     norm_tign_new = norms(k,1);
 308
 309
 310     %only do norm for times before final detection time
 311     time_mask = tign_new < max(tign_new(:)); %max(pts(end,3));  %max(max(pts(:,3)));
 312     norms(k,2) = norm(tign_new(time_mask)-tign_old(time_mask))/norm(tign_old(time_mask)-min(tign_old(:)));
 313     figure(fig_num+3);plot(norms(1:k,2));
 314     xlabel('Iterations of Interpolation')
 315     ylabel('Relative of difference')
 316     tstr = sprintf('Relative difference between successive \n TIGN after each interpolation');
 317     title(tstr)
 318     figure(fig_num+4);plot(1:k,norms(1:k,1))
 319     tstr = sprintf('Norm of difference between times of \n detections and TIGN at detection locations');
 320     title(tstr);
 321     xlabel('Iterations of Interpolation')
 322     ylabel('Norm of difference')
 323     %temp_var(k) = min(tign_new(:))-ps.red.start_datenum;
 324     %figure(fig_num+5);plot(1:k,temp_var(1:k)*24),title('Change in ignition time'),xlabel('iteration'),ylabel('hours')
 325     fprintf('Loop %d complete norm of diff = %f \n', k,norms(k))
 326     norm_break = 1;
 327     if k > 2 && norms(k,norm_break) >= norm(k-1,norm_break) && norms(k,norm_break) >= norms(k-2,norm_break)
 328         fprintf('graph norm increase \n')
 329         rm = 0;
 330         rt = 0;
 331         %break
 332     end
 333 end
 334 %tign_new = [];
 335 %smooth the result
 336
 337
 338 % %%% try using ground detections after paths....
 339 ground_steps = 6;
 340 tign_ground = tign_new;
 341 if gq
 342     data_area = sum(infire(:));
 343     for i = 1:ground_steps
 344         fire_mask = tign_ground < tmax-0.1;
 345         fire_area = sum(fire_mask(:));
 346         out_side = fire_mask(:)-infire(:);
 347         out_side(out_side<0)=0;
 348         out_sum(i) = sum(out_side(:));
 349         out_fraction = out_sum(i)/data_area;
 350         if i > 2 && out_sum(i) == out_sum(i-1)
 351             smooth_ground = 0.9*smooth_ground;
 352             fprintf('decreasing ground smoothing \n')
 353         end
 354         if out_fraction < 0.1
 355             %break
 356         end
 357         fprintf('Fire area: %f pixel_out/data area: %f  pixels_out %f\n',fire_area,out_sum(i)/data_area,out_sum(i));
 358         if use_beta_likes == 1
 359             tign_ground(~infire) = beta_vect(~infire).*tign_ground(~infire)+(1-beta_vect(~infire)).*tign_flat(~infire);
 360         else
 361             tign_ground(~infire) = beta*tign_ground(~infire)+(1-beta)*tign_flat(~infire);
 362         end
 363         %tign_ground(~infire) = beta*tign_ground(~infire)+(1-beta)*tign_flat(~infire);
 364         %tign_ground(~infire) = beta_vect(~infire).*tign_ground(~infire)+(1-beta_vect(~infire)).*tign_flat(~infire);
 365         t_mask = tign_ground > tmax;
 366         tign_ground(t_mask) = tmax;
 367         %tign_temp = imgaussfilt(tign_ground,1/6);
 368         %tign_temp = imgaussfilt(tign_ground,smooth_ground );
 369         tign_temp = smooth_up(tign_ground,a,b);
 370         tign_ground(~infire) = tign_temp(~infire);
 371 %         a = 1/10;%1/2-1/(2*i);
 372 %         tign_ground(infire) = a*tign_ground(infire)+(1-a)*tign_temp(infire);
 373 %         figure(73),scatter3(ps.red.fxlong(~infire),ps.red.fxlat(~infire),tign_ground(~infire))
 374 %         pause(1/2)
 375 % plot perimeter shrinking
 376 %         figure(159)
 377 %         contourf(ps.red.fxlong,ps.red.fxlat,tign_ground,20,'k'),hold on
 378 %         scatter(pts(1:5:end,2),pts(1:5:end,1),'*r'),hold off
 379 %         t_str = sprintf('Perimeter Shrinking \n Iteration %d',i);
 380 %         save_str = sprintf('perim_shrink_%d',i);
 381 %         xlabel('Lon [degree]'),ylabel('Lat [degre]'),title(t_str)
 382 %         close 159
 383 %         savefig(save_str);
 384 %         saveas(gcf,[save_str '.png']);
 385         %figure(160),mesh(ps.red.fxlong,ps.red.fxlat,tign_ground);
 386         %pause(.5)
 387         %t_min(i) = min(tign_ground(:));
 388     end
 389 end
 390 tign_new = tign_ground;
 391 %% end of post smoothing
 392
 393 tign_new = smooth_up(tign_new,a,b);
 394
 395
 396
 397 % %fix the ripple on the top induced by the smoothing
 398 r_mask = tign_new >=max(tign_new(:))-0.05;
 399 tign_new(r_mask) = max(tign_new(:));
 400
 401
 402 figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new-t0)
 403 title(title_str)
 404 hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3)-t0,'*r'),hold off
 405
 406 %plot the forecast for comparison
 407 figure,mesh(ps.red.fxlong,ps.red.fxlat,ps.red.tign-t0)
 408 title('Forecast')
 409 hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3)-t0,'*r')
 410
 411 filled_contours(ps,tign_new)
 412
 413 end
 414
 415