detection/detect_fit_linesearch.m

   1 function p=detect_fit_linesearch(prefix)
   2 % from a copy of barker2
   3
   4 disp('input data')
   5     % to create conus.kml:
   6     % download http://firemapper.sc.egov.usda.gov/data_viirs/kml/conus_hist/2012/conus_20120914.kmz
   7     % and gunzip
   8     %
   9     % to create w.mat:
  10     % run Adam's simulation, currently results in
  11     % /share_home/akochans/NASA_WSU/wrf-fire/WRFV3/test/em_barker_moist/wrfoutputfiles_live_0.25
  12     % then in Matlab
  13     % f='wrfout_d05_2012-09-15_00:00:00';
  14     % t=nc2struct(f,{'Times'},{});  n=size(t.times,2);
  15     % w=nc2struct(f,{'TIGN_G','FXLONG','FXLAT','UNIT_FXLAT','UNIT_FXLONG','Times',NFUEL_CAT'},{'DX','DY'},n);
  16     % save ~/w.mat w
  17     %
  18     % fuels.m is created by WRF-SFIRE at the beginning of the run
  19
  20     % ****** REQUIRES Matlab 2013a - will not run in earlier versions *******
  21
  22
  23     v=read_fire_kml('conus_viirs.kml');
  24     detection='VIIRS';
  25     if ~exist('prefix','var'),
  26         prefix='viirs';
  27     end
  28
  29     a=load('w');w=a.w;
  30     if ~isfield(w,'dx'),
  31         w.dx=444.44;
  32         w.dy=444.44;
  33         warning('fixing up w for old w.mat file from Barker fire')
  34     end
  35     dx=w.dx;
  36     dy=w.dy;
  37
  38     fuel.weight=0; % just to let Matlab know what fuel is going to be at compile time
  39     fuels
  40
  41
  42 disp('subset and process inputs')
  43
  44     % establish boundaries from simulations
  45
  46     sim.min_lat = min(w.fxlat(:));
  47     sim.max_lat = max(w.fxlat(:));
  48     sim.min_lon = min(w.fxlong(:));
  49     sim.max_lon = max(w.fxlong(:));
  50     sim.min_tign= min(w.tign_g(:));
  51     sim.max_tign= max(w.tign_g(:));
  52
  53     % active
  54     act.x=find(w.tign_g(:)<sim.max_tign);
  55     act.min_lat = min(w.fxlat(act.x));
  56     act.max_lat = max(w.fxlat(act.x));
  57     act.min_lon = min(w.fxlong(act.x));
  58     act.max_lon = max(w.fxlong(act.x));
  59
  60     % domain bounds
  61     margin=0.3;
  62     fprintf('enter relative margin around the fire (%g)',margin);
  63     in=input(' > ');
  64     if ~isempty(in),margin=in;end
  65     dis.min_lon=max(sim.min_lon,act.min_lon-margin*(act.max_lon-act.min_lon));
  66     dis.min_lat=max(sim.min_lat,act.min_lat-margin*(act.max_lat-act.min_lat));
  67     dis.max_lon=min(sim.max_lon,act.max_lon+margin*(act.max_lon-act.min_lon));
  68     dis.max_lat=min(sim.max_lat,act.max_lat+margin*(act.max_lat-act.min_lat));
  69
  70     default_bounds{1}=[sim.min_lon,sim.max_lon,sim.min_lat,sim.max_lat];
  71     descr{1}='fire domain';
  72     default_bounds{2}=[dis.min_lon,dis.max_lon,dis.min_lat,dis.max_lat];
  73     descr{2}='around fire';
  74     default_bounds{3}=[-119.5, -119.0, 47.95, 48.15];
  75     descr{3}='Barker fire';
  76     for i=1:length(default_bounds),
  77         fprintf('%i: %s %8.5f %8.5f %8.5f %8.5f\n',i,descr{i},default_bounds{i});
  78     end
  79     bounds=input_num('bounds [min_lon,max_lon,min_lat,max_lat] or number of bounds above',2);
  80     if length(bounds)==1,
  81         bounds=default_bounds{bounds};
  82     end
  83     fprintf('using bounds %8.5f %8.5f %8.5f %8.5f\n',bounds)
  84     display_bounds=bounds;
  85
  86     [ii,jj]=find(w.fxlong>=bounds(1) & w.fxlong<=bounds(2) & w.fxlat >=bounds(3) & w.fxlat <=bounds(4));
  87     ispan=min(ii):max(ii);
  88     jspan=min(jj):max(jj);
  89
  90     % restrict data
  91     fxlat=w.fxlat(ispan,jspan);
  92     fxlong=w.fxlong(ispan,jspan);
  93     tign_g=w.tign_g(ispan,jspan);
  94     nfuel_cat=w.nfuel_cat(ispan,jspan);
  95
  96     min_lon = display_bounds(1);
  97     max_lon = display_bounds(2);
  98     min_lat = display_bounds(3);
  99     max_lat = display_bounds(4);
 100
 101     % convert tign_g to datenum as tign, based zero at the end
 102     % assuming there is some place not on fire yet where tign_g = w.times
 103     %
 104     w_time_datenum=datenum(char(w.times)'); % the timestep of the wrfout, in days
 105     max_sim_time=max(tign_g(:));          % max time in the simulation, in sec
 106     tign=(tign_g - max_sim_time)/(24*60*60) + w_time_datenum; % assume same
 107
 108     % tign_g = max_sim_time + (24*60*60)*(tign - w_time_datenum)
 109     min_tign=min(tign(:));
 110     max_tign=max(tign(:));
 111
 112     % rebase time on the largest tign_g = the time of the first frame with fire, in days
 113     base_time=min_tign;
 114
 115     v.tim = v.tim - base_time;
 116     tign= tign - base_time;
 117
 118     % select fire detection within the domain and time
 119     bii=(v.lon > min_lon & v.lon < max_lon & v.lat > min_lat & v.lat < max_lat);
 120
 121     tol=0.01;
 122     tim_in = v.tim(bii);
 123     u_in = unique(tim_in);
 124     fprintf('detection times from ignition\n')
 125     for i=1:length(u_in)
 126         detection_freq(i)=sum(tim_in>u_in(i)-tol & tim_in<u_in(i)+tol);
 127         fprintf('%8.5f days %s UTC %3i %s detections\n',u_in(i),...
 128         datestr(u_in(i)+base_time),detection_freq(i),detection);
 129     end
 130     [max_freq,i]=max(detection_freq);
 131 %    detection_bounds=input_num('detection bounds as [upper,lower]',...
 132 %        [u_in(i)-min_tign-tol,u_in(i)-min_tign+tol]);
 133     detection_bounds = [u_in(i)-tol,u_in(i)+tol];
 134     bi = bii & detection_bounds(1) <= v.tim & v.tim <= detection_bounds(2);
 135     % now detection selected in time and space
 136     lon=v.lon(bi);
 137     lat=v.lat(bi);
 138     res=v.res(bi);
 139     tim=v.tim(bi);
 140     tim_ref = mean(tim);
 141
 142     fprintf('%i detections selected\n',sum(bi))
 143     detection_time=tim_ref;
 144     detection_datenum=tim_ref+base_time;
 145     detection_datestr=datestr(tim_ref+base_time);
 146     fprintf('mean detection time %g days from ignition %s UTC\n',...
 147         detection_time,detection_datestr);
 148     fprintf('days from ignition  min %8.5f max %8.5f\n',min(tim)-min_tign,max(tim)-min_tign);
 149     fprintf('longitude           min %8.5f max %8.5f\n',min(lon),max(lon));
 150     fprintf('latitude            min %8.5f max %8.5f\n',min(lat),max(lat));
 151
 152     % set up reduced resolution plots
 153     [m,n]=size(fxlong);
 154     m_plot=m; n_plot=n;
 155
 156     m1=map_index(display_bounds(1),bounds(1),bounds(2),m);
 157     m2=map_index(display_bounds(2),bounds(1),bounds(2),m);
 158     n1=map_index(display_bounds(3),bounds(3),bounds(4),n);
 159     n2=map_index(display_bounds(4),bounds(3),bounds(4),n);
 160     mi=m1:ceil((m2-m1+1)/m_plot):m2; % reduced index vectors
 161     ni=n1:ceil((n2-n1+1)/n_plot):n2;
 162     mesh_fxlong=fxlong(mi,ni);
 163     mesh_fxlat=fxlat(mi,ni);
 164     [mesh_m,mesh_n]=size(mesh_fxlat);
 165
 166     % find ignition point
 167     [i_ign,j_ign]=find(tign == min(tign(:)));
 168     if length(i_ign)~=1,error('assuming single ignition point here'),end
 169
 170     % set up constraint on ignition point being the same
 171     Constr_ign = zeros(m,n); Constr_ign(i_ign,j_ign)=1;
 172
 173     %
 174     % *** create detection mask for data likelihood ***
 175     %
 176     detection_mask=zeros(m,n);
 177     detection_time=tim_ref*ones(m,n);
 178
 179     % resolution diameter in longitude/latitude units
 180     rlon=0.5*res/w.unit_fxlong;
 181     rlat=0.5*res/w.unit_fxlat;
 182
 183
 184     lon1=lon-rlon;
 185     lon2=lon+rlon;
 186     lat1=lat-rlat;
 187     lat2=lat+rlat;
 188     for i=1:length(lon),
 189         square = fxlong>=lon1(i) & fxlong<=lon2(i) & ...
 190                  fxlat >=lat1(i) & fxlat <=lat2(i);
 191         detection_mask(square)=1;
 192     end
 193
 194     % for display in plotstate
 195     C=0.5*ones(1,length(res));
 196     X=[lon1,lon2,lon2,lon1]';
 197     Y=[lat1,lat1,lat2,lat2]';
 198 %    plotstate(1,detection_mask,['Fire detection at ',detection_datestr],[])
 199     % add ignition point
 200 %    hold on, plot(w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign),'xw'); hold off
 201     % legend('first ignition at %g %g',w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign))
 202
 203     fuelweight(length(fuel)+1:max(nfuel_cat(:)))=NaN;
 204     for j=1:length(fuel),
 205         fuelweight(j)=fuel(j).weight;
 206     end
 207     W = zeros(m,n);
 208     for j=1:n, for i=1:m
 209            W(i,j)=fuelweight(nfuel_cat(i,j));
 210     end,end
 211
 212 %    plotstate(2,W,'Fuel weight',[])
 213
 214 disp('optimization loop')
 215 h =zeros(m,n); % initial increment
 216 plotstate(3,tign,'Forecast fire arrival time',detection_time(1));
 217 print('-dpng','tign_forecast.png');
 218
 219 forecast=tign;
 220 mesh_tign_detect(4,fxlong,fxlat,forecast,v,'Forecast fire arrival time')
 221
 222 fprintf('********** Starting iterations **************\n');
 223
 224 % can change the objective function here
 225 alpha=input_num('penalty coefficient alpha',1/1000);
 226 if(alpha < 0)
 227     error('Alpha is not allowed to be negative.')
 228 end
 229
 230 % TC = W/(900*24); % time constant = fuel gone in one hour
 231 TC = 1/24;  % detection time constants in hours
 232 stretch=input_num('Tmin,Tmax,Tneg,Tpos',[0.5,10,5,10]);
 233 nodetw=input_num('no fire detection weight',0.5);
 234 power=input_num('negative laplacian power',1.02);
 235
 236 % storage for h maps
 237 maxiter = 2;
 238 maxdepth=2;
 239 h_stor = zeros(m,n,maxiter);
 240
 241 for istep=1:maxiter
 242
 243     fprintf('********** Iteration %g/%g **************\n', istep, maxiter);
 244
 245     psi = detection_mask - nodetw*(1-detection_mask);
 246
 247     % initial search direction, normed so that max(abs(search(:))) = 1.0
 248     [Js,search]=objective(tign,h);
 249     search = -search/big(search);
 250
 251     plotstate(5,search,'Search direction',0);
 252     print('-dpng', sprintf('%s_search_dir_%d.png', prefix, istep));
 253     [Jsbest,best_stepsize] = linesearch(4.0,Js,tign,h,search,4,maxdepth);
 254 %    plotstate(21,tign+h+3*search,'Line search (magic step_size=3)',detection_time(1));
 255     fprintf('Iteration %d: best step size %g\n', istep, best_stepsize);
 256     if(best_stepsize == 0)
 257         disp('Cannot improve in this search direction anymore, exiting now.');
 258         break;
 259     end
 260     h = h + best_stepsize*search;
 261     plotstate(10+istep,tign+h,sprintf('Analysis iteration %i [Js=%g]',istep,Jsbest),detection_time(1));
 262     print('-dpng',sprintf('%s_descent_iter_%d.png', prefix, istep));
 263     h_stor(:,:,istep) = h;
 264 end
 265 % rebase the analysis to the original simulation time
 266 analysis=tign+h;
 267 % w.tign_g = max_sim_time + (24*60*60)*(tign - w_time_datenum)
 268
 269 mesh_tign_detect(6,fxlong,fxlat,analysis,v,'Analysis fire arrival time')
 270 mesh_tign_detect(7,fxlong,fxlat,analysis-forecast,[],'Analysis - forecast difference')
 271
 272 [p.red.tign,p.red.tign_datenum] = rebase_time_back(tign+h);
 273 % analysis = max_sim_time + (24*60*60)*(tign+h + base_time - w_time_datenum);
 274 % err=big(p.tign_sim-analysis)
 275 [p.time.sfire,p.time.datenum] = rebase_time_back(detection_bounds);
 276 p.time.datestr=datestr(p.time.datenum);
 277 p.tign_g=w.tign_g;
 278 p.tign_g(ispan,jspan)=p.red.tign;
 279
 280 % max_sim_time + (24*60*60)*(tign+h + base_time - w_time_datenum);
 281
 282 disp('input the analysis as tign in WRF-SFIRE with fire_perimeter_time=detection time')
 283
 284 figure(9);
 285 col = 'rgbck';
 286 fill(X,Y,C,'EdgeAlpha',1,'FaceAlpha',0);
 287 for j=1:maxiter
 288     contour(mesh_fxlong,mesh_fxlat,tign+h_stor(:,:,j),[detection_time(1),detection_time(1)],['-',col(j)]); hold on
 289 end
 290 hold off
 291 title('Contour changes vs. step');
 292 xlabel('Longitude');
 293 ylabel('Latitude');
 294 print('-dpng',sprintf( '%s_contours.png', prefix));
 295
 296     function [time_sim,time_datenum]=rebase_time_back(time_in)
 297         time_datenum = time_in + base_time;
 298         time_sim = max_sim_time + (24*60*60)*(time_datenum - w_time_datenum);
 299     end
 300
 301     function varargout=objective(tign,h,doplot)
 302         % [J,delta]=objective(tign,h,doplot)
 303         % J=objective(tign,h,doplot)
 304         % compute objective function and optionally gradient delta direction
 305         T=tign+h;
 306         [f0,f1]=like1(psi,detection_time-T,TC*stretch);
 307         F = f1;             % forcing
 308         % objective function and preconditioned gradient
 309         Ah = poisson_fft2(h,[dx,dy],power);
 310         % compute both parts of the objective function and compare
 311         J1 = 0.5*(h(:)'*Ah(:));
 312         J2 = -ssum(f0);
 313         J = alpha*J1 + J2;
 314         fprintf('Objective function J=%g (J1=%g, J2=%g)\n',J,J1,J2);
 315         if nargout==1,
 316             varargout={J};
 317             return
 318         end
 319         gradJ = alpha*Ah + F;
 320         fprintf('Gradient: norm Ah %g norm F %g\n', norm(Ah,2), norm(F,2));
 321         if exist('doplot','var'),
 322             plotstate(7,f0,'Detection likelihood',0.5,'-w');
 323             plotstate(8,f1,'Detection likelihood derivative',0);
 324             plotstate(9,F,'Forcing',0);
 325             plotstate(10,gradJ,'gradient of J',0);
 326         end
 327         delta = solve_saddle(Constr_ign,h,F,0,@(u) poisson_fft2(u,[dx,dy],-power)/alpha);
 328         varargout=[{J},{delta}];
 329         % plotstate(11,delta,'Preconditioned gradient',0);
 330         %fprintf('norm(grad(J))=%g norm(delta)=%g\n',norm(gradJ,'fro'),norm(delta,'fro'))
 331     end
 332
 333     function plotstate(fig,T,s,c,linespec)
 334         fprintf('Figure %i %s\n',fig,s)
 335         plotmap(fig,mesh_fxlong,mesh_fxlat,T(mi,ni),' ');
 336         hold on
 337         hh=fill(X,Y,C,'EdgeAlpha',1,'FaceAlpha',0);
 338         if ~exist('c','var') || isempty(c) || isnan(c),
 339             title(s);
 340         else
 341             title(sprintf('%s, contour=%g',s,c(1)))
 342             if ~exist('linespec','var'),
 343                 linespec='-k';
 344             end
 345             contour(mesh_fxlong,mesh_fxlat,T(mi,ni),[c c],linespec)
 346         end
 347         hold off
 348         ratio=[w.unit_fxlat,w.unit_fxlong];
 349         xlabel longtitude
 350         ylabel latitude
 351         ratio=[ratio/norm(ratio),1];
 352         daspect(ratio)
 353         axis tight
 354         drawnow
 355     end
 356
 357
 358     function [Jsmin,best_stepsize] = linesearch(max_step,Js0,tign,h,search,nmesh,max_depth)
 359         step_low = 0;
 360         Jslow = Js0;
 361         step_high = max_step;
 362         Jshigh = objective(tign,h+max_step*search);
 363         for d=1:max_depth
 364             step_sizes = linspace(step_low,step_high,nmesh+2);
 365             Jsls = zeros(nmesh+2,1);
 366             Jsls(1) = Jslow;
 367             Jsls(nmesh+2) = Jshigh;
 368             for i=2:nmesh+1
 369                 Jsls(i) = objective(tign,h+step_sizes(i)*search);
 370             end
 371
 372             figure(8);
 373             plot(step_sizes,Jsls,'+-');
 374             title(sprintf('Objective function Js vs. step size, iter=%d,depth=%d',istep,d), 'fontsize', 16);
 375             xlabel('step\_size [-]','fontsize',14);
 376             ylabel('Js [-]','fontsize',14);
 377             print('-dpng',sprintf('%s_linesearch_iter_%d_depth_%d.png',prefix,istep,d));
 378
 379             [Jsmin,ndx] = min(Jsls);
 380
 381             low = max(ndx-1,1);
 382             high = min(ndx+1,nmesh+2);
 383             Jslow = Jsls(low);
 384             Jshigh = Jsls(high);
 385             step_low = step_sizes(low);
 386             step_high = step_sizes(high);
 387         end
 388
 389         best_stepsize = step_sizes(ndx);
 390     end
 391
 392 end % detect_fit
 393
 394 function i=map_index(x,a,b,n)
 395 % find image of x under linear map [a,b] -> [1,m]
 396 % and round to integer
 397 i=round(1+(n-1)*(x-a)/(b-a));
 398 end