cycling/squish.m

   1 function tign_new = squish(ps,ai)
   2 %% ps is struct with paths, red, graph,distances, etc.
   3 %% ps = graph_dets(w)
   4 %% fs - analysis or ignitions start, ai = 1 ==> analysis
   5 %             if ai == 1
   6 %                 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);
   7 %             end
   8
   9 %combine l2 time with detection points fixed to grid
  10 pts = ps.grid_pts;
  11 pts(:,3) = ps.points(:,3);
  12 tign = ps.red.tign;
  13 [m,n] = size(tign);
  14 idx = ps.idx;
  15 fig_num = 23;
  16 pts_length = length(ps.grid_pts);
  17 max_l2_time = max(max(pts(:,3)));
  18 if ai == 1
  19     fig_num = fig_num+1;
  20     tign_new = ps.red.tign;
  21     tign_flat=ps.red.end_datenum*ones(size(tign));
  22     title_str = 'Analysis';
  23 else
  24     %max_l2_time = max(tign(:));
  25     %flat tign_new
  26     f_days=input_num('Days of forecast to produce? ',0);
  27     tign_new=max_l2_time*ones(size(tign))+f_days;
  28     if f_days > 0
  29         ps.red.end_datenum = ps.red.end_datenum + f_days;
  30     end
  31     %tign_new = estimate_tign(ps);
  32     figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new)
  33     title_str = 'Generated fire cone';
  34 end
  35
  36
  37 %mask = ones(size(tign));
  38 %change back to z = tign
  39 % if exist('ps.tn','var')
  40 %     figure(68),contour3(ps.red.fxlong,ps.red.fxlat,ps.tn-ps.red.start_datenum,20)
  41 % else
  42 %     figure(68),contour3(ps.red.fxlong,ps.red.fxlat,tign-ps.red.start_datenum,20)
  43 % end
  44 % hold on
  45 %lons = ps.red.fxlong;
  46 %lats = ps.red.fxlat;
  47
  48 %make matrix of tign time for the points in the graph
  49 t_times = zeros(pts_length,1);
  50 for i = 1:pts_length
  51     t_times(i) = ps.red.tign(ps.idx(i,1),ps.idx(i,2));
  52 end
  53
  54 %just squish the ignition point
  55 %bigger k ==> more smoothing
  56
  57 norms=[];
  58
  59 %random multiplier, increase for larger grids
  60 %perturbs points on path in x-y plane
  61 % try computing this as a fraction of grid size
  62
  63 rm = 2;
  64
  65 % random multiplier, keep the same
  66 % perturbs points downward in time to
  67 rt = 0.5;
  68 % weight for tign_new
  69
  70 %alhpa blends  estimate of tign at a point with old estimate
  71 % new_estimate = alph*current_setimate + (1-alpha)*old_estimate
  72 alpha = 0.5;
  73 %constant for smooth in rls_shp
  74 alpha_2 = 0.05; %smaller alph_2 ==> smoother
  75 %number of loops to run
  76 smoothings = 40;
  77 for k = 1:smoothings
  78     rm = max(round(smoothings/k),2);
  79     figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new)
  80     title(title_str)
  81     hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3),'*r'),hold off
  82     pause(1/k)
  83     for i = 1:length(ps.paths)
  84         p = ps.paths(i).p;
  85         %         figure(73),hold on
  86         %         plot3(pts(p,2),pts(p,1),pts(p,3)-ps.red.start_datenum,'r')
  87         %         hold off
  88         %plot3(pts(p,2),pts(p,1),tign(idx(p,1),idx(p,2))-ps.red.start_datenum,'g')
  89         for j = 1:length(p)
  90             tign_old = tign_new;
  91             %mesh indices for path points, perturbed
  92             p_i = idx(p(j),1);%+rm*round(randn);
  93             p_j = idx(p(j),2);%+rm*round(randn);
  94             near_bound = p_i-rm > 0 && p_i + rm < m & p_j + rm > 0 && p_j < n;
  95             %%% make mean of old and new, in small ransdom size block around path point
  96             if ~near_bound
  97                 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*rand;
  98             else %block of a signle pointmake small
  99                 tign_new(p_i,p_j) = alpha*tign_new(p_i,p_j) + (1-alpha)*pts(p(j),3)-rt*rand;
 100             end
 101             %%%% alternate strategy., not working
 102             %             if k == 1 && j > 1
 103             %             t1 = min(tign_new(p_i,p_j),tign_new(p_i-1,p_j-1));
 104             %             t2 = max(tign_new(p_i,p_j),tign_new(p_i-1,p_j-1));
 105             %             dt12 = t2-t1;
 106             %             dg = pts(p(j),3)-pts(p(j-1),3);
 107             %             time_shift = 0.5*(dg-dt12);
 108             %             t1 = t1-time_shift;
 109             %             t2 = t2+time_shift;
 110             %             tign_new(p_i,p_j) = max(t1,t2);
 111             %             tign_new(p_i-1,p_j-1) = min(t1,t2);
 112             %             end
 113             if ai == 1
 114                 if ~near_bound
 115                     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);
 116                 else
 117                     tign_flat(p_i,p_j) = 0.5*(tign_new(p_i,p_j) + pts(p(j),3)-rt*rand);
 118                 end
 119             end
 120
 121             % interpolate new point between adjacent path detections
 122             if j > 1
 123                 %weighted average will move new point close to that with
 124                 %higher FRP
 125                 %                 frp1 = ps.points(p(j-1),5);
 126                 %                 frp2 = ps.points(p(j),5);
 127                 %                 w1 = frp1/(frp1+frp2);
 128                 %                 w2 = frp2/(frp1+frp2);
 129                 w1 = 0.5;
 130                 w2 = 0.5;
 131                 %assign tign for all in small, block around midpoint
 132                 %             frp1 = ps.points(p(j-1),5);
 133                 %             frp2 = ps.points(p(j),5);
 134                 %             w1 = frp1/(frp1+frp2);
 135                 %             w2 = frp2/(frp1+frp2);
 136                 new_lat = w1*pts(p(j-1),1)+w2*pts(p(j),1);
 137                 new_lon = w1*pts(p(j-1),2)+w2*pts(p(j),2);
 138                 new_t = w1*pts(p(j-1),3)+w2*pts(p(j),3);
 139                 [new_i,new_j,new_lat,new_lon]= fixpt(ps.red,[new_lat,new_lon]);
 140                 %flatten area around new pt
 141                 %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;
 142                 % why do it again?
 143                 if ai == 1
 144                     if ~near_bound
 145                         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;
 146                     end
 147                 end
 148                 %mask(new_i,new_j)=0;
 149             end
 150         end
 151     end
 152     %size of local averaging to apply aoutomate by grid size?
 153     patch = 6;
 154
 155     %smooth the tign
 156     %tign_new = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_new);
 157     tign_new = rlx_shp(tign_new,alpha_2,patch);
 158     if ai == 1
 159         tign_flat = rlx_shp(tign_flat,alpha_2,patch);
 160         %tign_flat = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_flat);
 161     end
 162     %collect information about tign at the detection points
 163     for i = 1:pts_length
 164         t_times(i) = tign_new(ps.idx(i,1),ps.idx(i,2));
 165         if ai ==1
 166             flat_times(i) = tign_flat(ps.idx(i,1),ps.idx(i,2));
 167         end
 168     end
 169     norms(k,1) = norm(t_times-ps.points(:,3),2);
 170     % blend flat start with forecast start for analysis
 171     if ai == 1
 172         %tign_flat = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_flat);
 173         tign_flat = rlx_shp(tign_flat,alpha_2,patch);
 174         norm_flat = norm(flat_times-ps.points(:,3),2);
 175         norm_tign_new = norms(k,1);
 176         if (norm_flat > norm_tign_new)
 177             %fprintf('blending flat start with forecast \n')
 178             tign_new = 0.5*(tign_new+tign_flat);
 179         end
 180     end
 181     %only do norm for times before final detection time
 182     time_mask = tign_new < pts(end,3);  %max(max(pts(:,3)));
 183     norms(k,2) = norm(tign_new(time_mask)-tign_old(time_mask));
 184 %     figure(fig_num+3);plot(norms(1:k,1));
 185 %     tstr = sprintf('Norm of difference between successive \n TIGN after each interpolation');
 186 %     title(tstr)
 187 %     figure(fig_num+4);plot(1:k,norms(1:k,2))
 188 %     tstr = sprintf('Norm of difference between times of \n detections and TIGN at detectionon locations');
 189 %     title(tstr);
 190 %     xlabel('Iterations of Interpolation')
 191     temp_var(k) = min(tign_new(:))-ps.red.start_datenum;
 192     %figure(fig_num+5);plot(1:k,temp_var(1:k)*24),title('Change in ignition time'),xlabel('iteration'),ylabel('hours')
 193     fprintf('Loop %d complete norm of diff = %f \n', k,norms(k))
 194     if k > 2 && norms(k,1) > norm(k-1,1) && norms(k,1) > norms(k-2,1)
 195         fprintf('graph norm increase \n')
 196         break
 197     end
 198 end
 199 %final smoothin after iterations
 200 tign_new = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_new);
 201 figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new)
 202     title(title_str)
 203     hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3),'*r'),hold off
 204 %tign_new = [];
 205 if ai == 1
 206     figure(fig_num+2),mesh(ps.red.fxlong,ps.red.fxlat,ps.red.tign)
 207     title('Forecast')
 208     hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3),'*r')
 209 end
 210 end
 211
 212 function new_mask = expand_mask(mask,border)
 213     [n,m] = size(mask);
 214     [ix,jx] = find(mask(border:n-border,border:m-border)==0);
 215     for i = 1:length(ix)
 216         mask(ix(i)-border:ix(i)+border,jx(i)-border:jx(i)+border)=1-1/border;
 217     end
 218     new_mask = mask;
 219 end
 220
 221