detection/detect_fit_level2.m

   1 function a=detect_fit_level2(varargin)
   2 % a=detect_fit_level2
   3 % a=detect_fit_level2(cycle,time_bounds,disp_bounds,w,force)
   4 % arguments
   5 %   cycle       cycle number, for display and output file names
   6 %   time_bounds [observation start, end, spinup start, end] (datenum)
   7 %   disp_bounds [longitude min, max, latitude min,max] (degrees)
   8 %   w           the data structure if not to read from w.mat
   9 %   force       use defaults
  10
  11 % arguments
  12 cycle=[];       if nargin>=1,cycle=varargin{1};end
  13 time_bounds=[]; if nargin>=2,time_bounds=varargin{2};end
  14 disp_bounds=[]; if nargin>=3,disp_bounds=varargin{3};end
  15 w=[];           if nargin>=4,w=varargin{4};end
  16 force=0;;       if nargin>=5,force=varargin{5};end
  17 if nargin>5, error('too many arguments'),end
  18
  19 % to create w.mat:
  20
  21 % run Adam's simulation, currently results in
  22 % /share_home/akochans/WRF341F/wrf-fire/WRFV3/test/em_utfire_1d_med_4km_200m
  23 % then in Matlab:
  24 % set up paths by running setup.m in wrf-fire/WRFV3/test/em_fire
  25 % f='wrfout_d01_2013-08-20_00:00:00';
  26 % w=read_wrfout_tign(f);
  27 % save ~/w.mat w
  28 % fuels.m is created by WRF-SFIRE at the beginning of the run
  29 % copy w.mat and fuel.m to your machine where this fuction will run
  30 % start matlab, set up paths by running setup.m in wrf-fire/WRFV3/test/em_fire
  31 % converge active fires detection geotiff to matlab by
  32 % wrf-fire/other/Matlab/detection/geotiff2mat.py
  33
  34
  35 % ****** REQUIRES Matlab 2016a - will not run in earlier versions *******
  36
  37 % figures
  38 fig.fig_map=0;
  39 fig.fig_3d=0;
  40 fig.fig_interp=0;
  41
  42 plot_also_wind=0;
  43
  44 disp('Loading simulation')
  45
  46 if isempty(w),
  47     w=load('w');w=w.w;
  48 end
  49
  50 if plot_also_wind,
  51     % wind
  52     ss=load('ss');ss=ss.s;
  53     ss.steps=size(ss.times,2);
  54     ss.time=zeros(1,ss.steps);
  55     for i=1:ss.steps,
  56         ss.time(i)=datenum(char(ss.times(:,i)'));
  57     end
  58     ss.num=1:ss.steps;
  59     ss.min_time=min(ss.time);
  60     ss.max_time=max(ss.time);
  61     % interpolate surface wind to center of the grid
  62     ss.uh=0.5*(ss.uah(1:end-1,:,:)+ss.uah(2:end,:,:));
  63     ss.vh=0.5*(ss.vah(:,1:end-1,:)+ss.vah(:,2:end,:));
  64     fprintf('min wind field time   %s\n',datestr(ss.min_time));
  65     fprintf('max wind field time   %s\n',datestr(ss.max_time));
  66
  67 end
  68
  69 fuel=[]; fuels
  70
  71 disp('Subset simulation domain and convert time')
  72
  73 red=subset_domain(w,force);
  74 if ~isempty(disp_bounds)
  75     red.disp_bounds=disp_bounds;
  76 else
  77     red.disp_bounds=[red.min_lon,red.max_lon,red.min_lat,red.max_lat];
  78 end
  79 fprintf('display bounds %g %g %g %g\n',red.disp_bounds);
  80
  81 disp('Loading and subsetting detections')
  82 fire_choice = input_num('which fire? Patch: [0], Camp: [1], Cougar [3]',0,1)
  83 if fire_choice == 1
  84     prefix='../campTIFs/';
  85 elseif fire_choice == 0
  86     prefix='../TIFs/';
  87 else
  88     prefix = '../cougarTIFs/';
  89 end
  90 %prefix='../TIFpoint/';
  91 %prefix='../perimTIFs/';
  92 % the level2 file names processed by geotiff2mat.py
  93 p=sort_rsac_files(prefix);
  94
  95 if isempty(time_bounds),
  96     time_bounds=subset_detection_time(red,p);
  97 end
  98 print_time_bounds(red,'Simulation',red.start_datenum,red.end_datenum)
  99 print_time_bounds(red,'Detections',time_bounds(1),time_bounds(2))
 100 print_time_bounds(red,'Spinup    ',time_bounds(3),time_bounds(4))
 101 red.time_bounds=time_bounds;
 102
 103 gstr = sprintf('g_%d.mat',cycle);
 104 if ~exist(gstr,'file')
 105     g = load_subset_detections(prefix,p,red,time_bounds,fig);
 106     save(gstr,'g')
 107 else
 108     load(gstr)
 109     fprintf('Loaded existing g.mat file \n')
 110     if input_num('Enter [1] to load from scratch',1)
 111         g = load_subset_detections(prefix,p,red,time_bounds,fig);
 112         save(gstr,'g');
 113     end
 114 end
 115
 116 [m,n]=size(red.fxlong);
 117
 118 % find ignition point
 119 tign=red.tign;
 120 %%%% remove this
 121 %tign(478,383)=tign(478,383)-1e4;
 122 [i_ign,j_ign]=find(tign ==min(tign(:)));
 123 if length(i_ign)~=1,error('assuming single ignition point here'),end
 124
 125 % set up constraint on ignition point being the same
 126 params.Constr_ign = zeros(m,n); params.Constr_ign(i_ign,j_ign)=1;
 127
 128
 129 % Parameters of the objective function
 130 params.alpha=input_num('penalty coefficient alpha',1/1000,force);
 131 % TC = W/(900*24); % time constant = fuel gone in one hour
 132 params.TC = 1/24;  % detection time constants in hours
 133 params.stretch=input_num('Tmin,Tmax,Tneg,Tpos',[0.5,10,5,10],force);
 134 params.weight=input_num('water,land,low,nominal,high confidence fire',[-1,-1,0.2,0.6,1],force);
 135 params.power=input_num('correction smoothness',1.02,force);
 136 params.doplot=0;
 137 params.dx=444;
 138 params.dy=444;
 139
 140 forecast=tign;
 141
 142 disp('optimization loop')
 143 h =zeros(m,n); % initial increment
 144 plot_state(3,red,'Forecast',forecast,g,time_bounds(1:2));
 145 savefig('forecast',cycle)
 146 plot_state_2d(4,red,'Forecast',forecast,g,time_bounds(2));
 147 savefig('forecast2d',cycle)
 148
 149
 150 fprintf('********** Starting iterations **************\n');
 151
 152
 153 % storage for h maps
 154 maxiter =2;
 155 maxdepth=3;
 156 h_stor = zeros(m,n,maxiter);
 157
 158 new_like = input_num('Use new likelihood? No == 0 ',0,1);
 159 for istep=1:maxiter
 160
 161     fprintf('********** Iteration %g/%g **************\n', istep, maxiter);
 162
 163     if new_like
 164         fprintf('Using new likelihood function\n')
 165         [p_like_spline,p_deriv_spline,n_deriv_spline] = make_spline(100,1000);
 166         %load splines.mat
 167         [Js,search]=detection_objective(tign,h,g,params,red,p_like_spline,p_deriv_spline,n_deriv_spline);
 168     else
 169         [Js,search]=detection_objective(tign,h,g,params,red);
 170     end
 171
 172
 173
 174     % initial search direction, normed so that max(abs(search(:))) = 1.0
 175     %[Js,search]=detection_objective(tign,h,g,params,red);
 176     search = -search/big(search);
 177
 178     print('-dpng', sprintf('%s_search_dir_%d.png', prefix, istep));
 179     [Jsbest,best_stepsize] = linesearch(4.0,Js,tign,h,search,4,maxdepth);
 180     plot_state(21,red,'Line search',tign+h+3*search,g,time_bounds(1:3));
 181     fprintf('Iteration %d: best step size %g\n', istep, best_stepsize);
 182     if(best_stepsize == 0)
 183         disp('Cannot improve in this search direction anymore, exiting now.');
 184         break;
 185     end
 186     h = h + best_stepsize*search;
 187
 188     %%   %dealing with the bump here
 189     bump_remove = 0;
 190     temp_analysis = tign+h;
 191     if bump_remove
 192         %figure,mesh(temp_analysis)
 193         corner_time = min([temp_analysis(1,1),temp_analysis(1,end),...
 194             temp_analysis(end,1),temp_analysis(end,end)]);
 195         if max(temp_analysis(:)) > corner_time
 196             fprintf('Bump forming \n');
 197             %figure,mesh(temp_analysis);
 198             mask_h = temp_analysis >= corner_time;
 199             h(mask_h) = 0;
 200             %figure,mesh(h);
 201         end
 202     end %bump removal
 203
 204     plot_state(10+istep,red,sprintf('Analysis iteration %i [Js=%g]',istep,Jsbest),tign+h,g,time_bounds(1:2));
 205     print('-dpng',sprintf('%s_descent_iter_%d.png', prefix, istep));
 206     h_stor(:,:,istep) = h;
 207 end
 208
 209 analysis=tign+h;
 210 %figure,mesh(analysis);
 211
 212 %cutting the top off of the "lid"
 213 bump_remove = 0;
 214 if bump_remove
 215     corner_time = min([temp_analysis(1,1),temp_analysis(1,end),...
 216         temp_analysis(end,1),temp_analysis(end,end)]);
 217     mask = analysis >= corner_time;
 218     analysis(mask) = corner_time;
 219 end
 220
 221 % w.tign_g = max_sim_time + (24*60*60)*(tign - w_time_datenum)
 222
 223 plot_state(6,red,'Analysis',analysis,g,time_bounds(1:2))
 224 savefig('analysis',cycle)
 225 plot_state_2d(5,red,'Analysis',analysis,g,time_bounds(2));
 226 savefig('analysis2d',cycle)
 227 plot_state_2d(7,red,{'Forecast','Analysis'},{forecast,analysis},g,time_bounds(2));
 228 savefig('forecast_analysis2d',cycle)
 229
 230 % spinup - combine analysis and forecast in the strip between
 231 % forecast fire area at restart time and outside of analysis fire area at perimeter time
 232
 233 restart_time=time_bounds(3);
 234 perimeter_time=time_bounds(4);
 235 wf = max(forecast - restart_time,0); % 0 in forecast fire area at restart time, >0 outside
 236 wa = max(perimeter_time-analysis,0); % 0 outside of analysis fire area at perimeter time, >0 inside
 237
 238 % check if we have inclusion so the strip exist
 239 shrink=nnz(wa + wf==0);
 240 if shrink,
 241     fprintf('Cannot spin up, fire area shrinking in analysis at %g points\n',shrink)
 242     error('Try an earlier restart time');
 243 end
 244
 245 % map the weights so that 0 ->1, 0->1,
 246 vf=1-wf./(wf+wa);  % 1  in forecast fire area at restart time, ADDING UP TO 1
 247 va=1-wa./(wf+wa); %  1  outside of analysis fire area at restart time, ADDING UP TO 1
 248
 249 % combine the forecast and analysis
 250 spinup = vf.*forecast + va.*analysis;
 251
 252 %figure,mesh(spinup),title('spinup mesh')
 253
 254 plot_state(8,red,'Spinup',spinup,g,time_bounds(3:4))
 255 savefig('spinup',cycle)
 256 plot_state(9,red,'Forecast for spinup',forecast,g,time_bounds(3:4))
 257 plot_state(10,red,'Analysis for spinup',analysis,g,time_bounds(3:4))
 258
 259
 260 % convert all time from datenum to seconds since run and insert in full
 261 % field
 262
 263 % insert analysis on reduced domain to the whole thing
 264 a.forecast=w.tign_g;
 265 a.analysis=w.tign_g;
 266 a.analysis(red.ispan,red.jspan)=datenum2time(analysis,red);
 267 a.spinup=w.tign_g;
 268 a.spinup(red.ispan,red.jspan)=datenum2time(spinup,red);
 269
 270 % display bounds
 271 a.disp_bounds=red.disp_bounds;
 272
 273 % copy time bounds to output structure
 274
 275 % as datenum - native
 276 a.observations_start_datenum=time_bounds(1);
 277 a.observations_end_datenum=time_bounds(2);
 278 a.restart_datenum=time_bounds(3);
 279 a.fire_perimeter_datenum=time_bounds(4);
 280
 281 % as seconds since start
 282 a.observations_start_time=datenum2time(time_bounds(1),red);
 283 a.observations_end_time=datenum2time(time_bounds(2),red);
 284 a.restart_time=datenum2time(time_bounds(3),red);
 285 a.fire_perimeter_time=datenum2time(time_bounds(4),red);
 286
 287 % as date strings
 288 a.observations_start_Times=stime(time_bounds(1),red);
 289 a.observations_end_Times=stime(time_bounds(2),red);
 290 a.restart_Times=stime(time_bounds(3),red);
 291 a.fire_perimeter_Times=stime(time_bounds(4),red);
 292
 293 % as days since start
 294 a.observations_start_days=a.observations_start_time/(24*3600);
 295 a.observations_end_days=a.observations_end_time/(24*3600);
 296 a.restart_days=a.restart_time/(24*3600);
 297 a.fire_perimeter_days=a.fire_perimeter_time/(24*3600);
 298
 299 a.cycle=cycle;
 300
 301 fprintf('Input the spinup as TIGN_G and restart\nfrom %s with fire_perimeter_time=%g\n',...
 302     a.restart_Times,a.fire_perimeter_time)
 303
 304 return
 305
 306     function [Jsmin,best_stepsize] = linesearch(max_step,Js0,tign,h,search,nmesh,max_depth)
 307         step_low = 0;
 308         Jslow = Js0;
 309         step_high = max_step;
 310         % Jshigh = detection_objective(tign,h+max_step*search,red);
 311         for d=1:max_depth
 312             step_sizes = linspace(step_low,step_high,nmesh+2);
 313             Jsls = zeros(nmesh+2,1);
 314             Jsls(1) = Jslow;
 315             % Jsls(nmesh+2) = Jshigh;
 316             for i=2:nmesh+2
 317                 %% possibly change h or tign here
 318                 Jsls(i) = detection_objective(tign,h+step_sizes(i)*search,g,params,red);
 319             end
 320             Jshigh=Jsls(nmesh+2);
 321             for i=1:nmesh+2
 322                 str=sprintf('step=%g objective function=%g',step_sizes(i),Jsls(i));
 323                 plot_state(20+i,red,str,tign+h+step_sizes(i)*search,g,time_bounds(1:2));
 324             end
 325
 326             figure(8);
 327             plot(step_sizes,Jsls,'+-');
 328             title(sprintf('Objective function Js vs. step size, iter=%d,depth=%d',istep,d), 'fontsize', 16);
 329             xlabel('step\_size [-]');
 330             ylabel('Js [-]');
 331             % print('-dpng',sprintf('%s_linesearch_iter_%d_depth_%d.png',prefix,istep,d));
 332
 333             [Jsmin,ndx] = min(Jsls);
 334
 335             low = max(ndx-1,1);
 336             high = min(ndx+1,nmesh+2);
 337             Jslow = Jsls(low);
 338             Jshigh = Jsls(high);
 339             step_low = step_sizes(low);
 340             step_high = step_sizes(high);
 341             if high<nmesh+2,
 342                 step_high = step_sizes(high);
 343             else
 344                 step_high = step_sizes(high)*2;
 345             end
 346         end
 347
 348         best_stepsize = step_sizes(ndx);
 349     end
 350
 351
 352 end % detect_fit_level2
 353
 354 function i=map_index(x,a,b,n)
 355 % find image of x under linear map [a,b] -> [1,m]
 356 % and round to integer
 357 i=round(1+(n-1)*(x-a)/(b-a));
 358 end
 359
 360 function savefig(file,cycle)
 361     if isempty(cycle),
 362         filename=file;
 363     else
 364         filename=sprintf('%s_%i',file,cycle);
 365     end
 366     h=gcf;
 367     fprintf('Saving figure %i as %s\n',h.Number,filename)
 368     print('-dpng',[filename,'.png'])
 369     saveas(gcf,[filename,'.fig'],'fig')
 370 end