cycling/squish3.m

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