other/Matlab/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,w)
   4 if nargin>=1,
   5     cycle=varargin{1};
   6 else
   7     cycle=[];
   8 end
   9 if nargin>=2,
  10     time_bounds=varargin{2};
  11 else
  12     time_bounds=[];
  13 end
  14 if nargin>=3,
  15     w=varargin{3};
  16 else
  17     w=[];
  18 end
  19
  20 % to create w.mat:
  21
  22 % run Adam's simulation, currently results in
  23 % /share_home/akochans/WRF341F/wrf-fire/WRFV3/test/em_utfire_1d_med_4km_200m
  24 % then in Matlab:
  25 % set up paths by running setup.m in wrf-fire/WRFV3/test/em_fire
  26 % f='wrfout_d01_2013-08-20_00:00:00';
  27 % t=nc2struct(f,{'Times'},{});  n=size(t.times,2)
  28 % w=nc2struct(f,{'Times','TIGN_G','FXLONG','FXLAT','UNIT_FXLAT','UNIT_FXLONG','XLONG','XLAT','NFUEL_CAT','ITIMESTEP'},{'DX','DY'.'DT'},n);
  29 % save ~/w.mat w
  30 % fuels.m is created by WRF-SFIRE at the beginning of the run
  31 % copy w.mat and fuel.m to your machine where this fuction will run
  32 % start matlab, set up paths by running setup.m in wrf-fire/WRFV3/test/em_fire
  33 % converge active fires detection geotiff to matlab by
  34 % wrf-fire/other/Matlab/detection/geotiff2mat.py
  35
  36
  37 % ****** REQUIRES Matlab 2016a - will not run in earlier versions *******
  38
  39 % figures
  40 fig.fig_map=0;
  41 fig.fig_3d=0;
  42 fig.fig_interp=0;
  43
  44 plot_also_wind=0;
  45
  46 disp('Loading simulation')
  47
  48 if isempty(w),
  49     w=load('w');w=w.w;
  50 end
  51
  52 if plot_also_wind,
  53     % wind
  54     ss=load('ss');ss=ss.s;
  55     ss.steps=size(ss.times,2);
  56     ss.time=zeros(1,ss.steps);
  57     for i=1:ss.steps,
  58         ss.time(i)=datenum(char(ss.times(:,i)'));
  59     end
  60     ss.num=1:ss.steps;
  61     ss.min_time=min(ss.time);
  62     ss.max_time=max(ss.time);
  63     % interpolate surface wind to center of the grid
  64     ss.uh=0.5*(ss.uah(1:end-1,:,:)+ss.uah(2:end,:,:));
  65     ss.vh=0.5*(ss.vah(:,1:end-1,:)+ss.vah(:,2:end,:));
  66     fprintf('min wind field time   %s\n',datestr(ss.min_time));
  67     fprintf('max wind field time   %s\n',datestr(ss.max_time));
  68
  69 end
  70
  71 fuel=[]; fuels
  72
  73 disp('Subset simulation domain and convert time')
  74
  75 red=subset_domain(w);
  76
  77 disp('Loading and subsetting detections')
  78
  79 prefix='TIFs/';
  80 % the level2 file names processed by geotiff2mat.py
  81 p=sort_rsac_files(prefix);
  82
  83 if isempty(time_bounds),
  84     time_bounds=subset_detection_time(red,p);
  85 end
  86
  87 g = load_subset_detections(prefix,p,red,time_bounds,fig);
  88
  89 [m,n]=size(red.fxlong);
  90
  91 % find ignition point
  92 tign=red.tign;
  93 [i_ign,j_ign]=find(tign == min(tign(:)));
  94 if length(i_ign)~=1,error('assuming single ignition point here'),end
  95
  96 % set up constraint on ignition point being the same
  97 params.Constr_ign = zeros(m,n); params.Constr_ign(i_ign,j_ign)=1;
  98
  99
 100 % Parameters of the objective function
 101 params.alpha=input_num('penalty coefficient alpha',1/1000);
 102 % TC = W/(900*24); % time constant = fuel gone in one hour
 103 params.TC = 1/24;  % detection time constants in hours
 104 params.stretch=input_num('Tmin,Tmax,Tneg,Tpos',[0.5,10,5,10]);
 105 params.weight=input_num('water,land,low,nominal,high confidence fire',[-1,-1,0.2,0.6,1]);
 106 params.power=input_num('correction smoothness',1.02);
 107 params.doplot=0;
 108 params.dx=444;
 109 params.dy=444;
 110
 111
 112 disp('optimization loop')
 113 h =zeros(m,n); % initial increment
 114 plot_state(3,red,'Forecast',tign,g,time_bounds(1:2));
 115 savefig('forecast',cycle)
 116
 117 forecast=tign;
 118
 119 fprintf('********** Starting iterations **************\n');
 120
 121
 122 % storage for h maps
 123 maxiter =2;
 124 maxdepth=3;
 125 h_stor = zeros(m,n,maxiter);
 126
 127 for istep=1:maxiter
 128
 129     fprintf('********** Iteration %g/%g **************\n', istep, maxiter);
 130
 131     % initial search direction, normed so that max(abs(search(:))) = 1.0
 132     [Js,search]=detection_objective(tign,h,g,params);
 133     search = -search/big(search);
 134
 135     print('-dpng', sprintf('%s_search_dir_%d.png', prefix, istep));
 136     [Jsbest,best_stepsize] = linesearch(4.0,Js,tign,h,search,4,maxdepth);
 137     plot_state(21,red,'Line search',tign+h+3*search,g,time_bounds(1:3));
 138     fprintf('Iteration %d: best step size %g\n', istep, best_stepsize);
 139     if(best_stepsize == 0)
 140         disp('Cannot improve in this search direction anymore, exiting now.');
 141         break;
 142     end
 143     h = h + best_stepsize*search;
 144     plot_state(10+istep,red,sprintf('Analysis iteration %i [Js=%g]',istep,Jsbest),...
 145         tign+h,g,time_bounds(1:2));
 146     print('-dpng',sprintf('%s_descent_iter_%d.png', prefix, istep));
 147     h_stor(:,:,istep) = h;
 148 end
 149 % rebase the analysis to the original simulation time
 150 analysis=tign+h;
 151 % w.tign_g = max_sim_time + (24*60*60)*(tign - w_time_datenum)
 152
 153 plot_state(6,red,'Analysis',analysis,g,time_bounds(1:2))
 154 savefig('analysis',cycle)
 155
 156 % spinup - combine analysis and forecast in the strip between
 157 % forecast fire area at restart time and outside of analysis fire area at perimeter time
 158
 159 restart_time=time_bounds(3);
 160 perimeter_time=time_bounds(4);
 161 wf = max(forecast - restart_time,0); % 0 in forecast fire area at restart time, >0 outside
 162 wa = max(perimeter_time-analysis,0); % 0 outside of analysis fire area at perimeter time, >0 inside
 163
 164 % check if we have inclusion so the strip exist
 165 shrink=nnz(wa + wf==0);
 166 if shrink,
 167     fprintf('Cannot spin up, fire area shrinking in analysis at %g points\n',shrink)
 168     error('Try an earlier restart time');
 169 end
 170
 171 % map the weights so that 0 ->1, 0->1,
 172 vf=1-wf./(wf+wa);  % 1  in forecast fire area at restart time, ADDING UP TO 1
 173 va=1-wa./(wf+wa); %  1  outside of analysis fire area at restart time, ADDING UP TO 1
 174
 175 % combine the forecast and analysis
 176 spinup = vf.*forecast + va.*analysis;
 177
 178 plot_state(8,red,'Spinup',spinup,g,time_bounds(3:4))
 179 savefig('spinup',cycle)
 180 plot_state(9,red,'Forecast for spinup',forecast,g,time_bounds(3:4))
 181 plot_state(10,red,'Analysis for spinup',analysis,g,time_bounds(3:4))
 182
 183
 184 % convert all time from datenum to seconds since run and insert in full
 185 % field
 186
 187 % insert analysis on reduced domain to the whole thing
 188 a.forecast=w.tign_g;
 189 a.analysis=w.tign_g;
 190 a.analysis(red.ispan,red.jspan)=datenum2time(analysis,red);;
 191 a.spinup=w.tign_g;
 192 a.spinup(red.ispan,red.jspan)=datenum2time(spinup,red);;
 193
 194 % copy time bounds to output structure
 195
 196 % as datenum - native
 197 a.observations_start_datenum=time_bounds(1);
 198 a.observations_end_datenum=time_bounds(2);
 199 a.restart_datenum=time_bounds(3);
 200 a.fire_perimeter_datenum=time_bounds(4);
 201
 202 % as seconds since start
 203 a.observations_start_time=datenum2time(time_bounds(1),red);
 204 a.observations_end_time=datenum2time(time_bounds(2),red);
 205 a.restart_time=datenum2time(time_bounds(3),red);
 206 a.fire_perimeter_time=datenum2time(time_bounds(4),red);
 207
 208 % as date strings
 209 a.observations_start_Times=stime(time_bounds(1),red);
 210 a.observations_end_Times=stime(time_bounds(2),red);
 211 a.restart_Times=stime(time_bounds(3),red);
 212 a.fire_perimeter_Times=stime(time_bounds(4),red);
 213
 214 % as days since start
 215 a.observations_start_days=a.observations_start_time/(24*3600);
 216 a.observations_end_days=a.observations_end_time/(24*3600);
 217 a.restart_days=a.restart_time/(24*3600);
 218 a.fire_perimeter_days=a.fire_perimeter_time/(24*3600);
 219
 220 fprintf('Input the spinup as TIGN_G and restart\nfrom %s with fire_perimeter_time=%g\n',...
 221     a.restart_Times,a.fire_perimeter_time)
 222
 223 return
 224
 225     function [Jsmin,best_stepsize] = linesearch(max_step,Js0,tign,h,search,nmesh,max_depth)
 226         step_low = 0;
 227         Jslow = Js0;
 228         step_high = max_step;
 229         % Jshigh = detection_objective(tign,h+max_step*search);
 230         for d=1:max_depth
 231             step_sizes = linspace(step_low,step_high,nmesh+2);
 232             Jsls = zeros(nmesh+2,1);
 233             Jsls(1) = Jslow;
 234             % Jsls(nmesh+2) = Jshigh;
 235             for i=2:nmesh+2
 236                 Jsls(i) = detection_objective(tign,h+step_sizes(i)*search,g,params);
 237             end
 238             Jshigh=Jsls(nmesh+2);
 239             for i=1:nmesh+2
 240                 str=sprintf('step=%g objective function=%g',step_sizes(i),Jsls(i));
 241                 plot_state(20+i,red,str,tign+h+step_sizes(i)*search,g,time_bounds(1:2));
 242             end
 243
 244             figure(8);
 245             plot(step_sizes,Jsls,'+-');
 246             title(sprintf('Objective function Js vs. step size, iter=%d,depth=%d',istep,d), 'fontsize', 16);
 247             xlabel('step\_size [-]');
 248             ylabel('Js [-]');
 249             % print('-dpng',sprintf('%s_linesearch_iter_%d_depth_%d.png',prefix,istep,d));
 250
 251             [Jsmin,ndx] = min(Jsls);
 252
 253             low = max(ndx-1,1);
 254             high = min(ndx+1,nmesh+2);
 255             Jslow = Jsls(low);
 256             Jshigh = Jsls(high);
 257             step_low = step_sizes(low);
 258             step_high = step_sizes(high);
 259             if high<nmesh+2,
 260                 step_high = step_sizes(high);
 261             else
 262                 step_high = step_sizes(high)*2;
 263             end
 264         end
 265
 266         best_stepsize = step_sizes(ndx);
 267     end
 268
 269
 270 end % detect_fit_level2
 271
 272 function i=map_index(x,a,b,n)
 273 % find image of x under linear map [a,b] -> [1,m]
 274 % and round to integer
 275 i=round(1+(n-1)*(x-a)/(b-a));
 276 end
 277
 278 function savefig(file,cycle)
 279     if isempty(cycle),
 280         filename=file;
 281     else
 282         filename=sprintf('%s_%i',file,cycle);
 283     end
 284     print('-dpng',[filename,'.png'])
 285     saveas(gcf,[filename,'.fig'],'fig')
 286 end