From 586f61d84d3ed7890b51f3c5b2fc2d184fc3b131 Mon Sep 17 00:00:00 2001 From: Jan Mandel Date: Fri, 24 Jun 2016 21:02:02 -0600 Subject: [PATCH] adding some graphs for detection_log_likelihood and fixing sign --- other/Matlab/detection/detect_fit_level2.m | 9 +++-- other/Matlab/detection/detection_log_likelihood.m | 17 ++++++-- other/Matlab/detection/detection_objective.m | 47 +++++++++++++++-------- other/Matlab/detection/plot_loglike.m | 22 +++++++++++ 4 files changed, 72 insertions(+), 23 deletions(-) create mode 100644 other/Matlab/detection/plot_loglike.m diff --git a/other/Matlab/detection/detect_fit_level2.m b/other/Matlab/detection/detect_fit_level2.m index 22f164cc..c31ce65a 100644 --- a/other/Matlab/detection/detect_fit_level2.m +++ b/other/Matlab/detection/detect_fit_level2.m @@ -88,6 +88,7 @@ end print_time_bounds(red,'Simulation',red.start_datenum,red.end_datenum) print_time_bounds(red,'Detections',time_bounds(1),time_bounds(2)) print_time_bounds(red,'Spinup ',time_bounds(3),time_bounds(4)) +red.time_bounds=time_bounds; g = load_subset_detections(prefix,p,red,time_bounds,fig); @@ -107,7 +108,7 @@ params.alpha=input_num('penalty coefficient alpha',1/1000); % TC = W/(900*24); % time constant = fuel gone in one hour params.TC = 1/24; % detection time constants in hours params.stretch=input_num('Tmin,Tmax,Tneg,Tpos',[0.5,10,5,10]); -params.weight=input_num('water,land,low,nominal,high confidence fire',[-1,-1,0.2,0.6,1]); % -1 ??? +params.weight=input_num('water,land,low,nominal,high confidence fire',[-10,-10,0.2,0.6,1]); % -1 ??? params.power=input_num('correction smoothness',1.02); params.doplot=0; params.dx=444; @@ -136,7 +137,7 @@ for istep=1:maxiter fprintf('********** Iteration %g/%g **************\n', istep, maxiter); % initial search direction, normed so that max(abs(search(:))) = 1.0 - [Js,search]=detection_objective(tign,h,g,params); + [Js,search]=detection_objective(tign,h,g,params,red); search = -search/big(search); print('-dpng', sprintf('%s_search_dir_%d.png', prefix, istep)); @@ -242,14 +243,14 @@ return step_low = 0; Jslow = Js0; step_high = max_step; - % Jshigh = detection_objective(tign,h+max_step*search); + % Jshigh = detection_objective(tign,h+max_step*search,red); for d=1:max_depth step_sizes = linspace(step_low,step_high,nmesh+2); Jsls = zeros(nmesh+2,1); Jsls(1) = Jslow; % Jsls(nmesh+2) = Jshigh; for i=2:nmesh+2 - Jsls(i) = detection_objective(tign,h+step_sizes(i)*search,g,params); + Jsls(i) = detection_objective(tign,h+step_sizes(i)*search,g,params,red); end Jshigh=Jsls(nmesh+2); for i=1:nmesh+2 diff --git a/other/Matlab/detection/detection_log_likelihood.m b/other/Matlab/detection/detection_log_likelihood.m index 3e4a8dc0..8ffe7dcb 100644 --- a/other/Matlab/detection/detection_log_likelihood.m +++ b/other/Matlab/detection/detection_log_likelihood.m @@ -3,23 +3,25 @@ function J=detection_log_likelihood(prefix,w) % in % prefix to search for files, such as 'TIFs/' % w file with extracted wrfout data, see detect_fit_level2 +% % out % J log likelhood -% figures +% figures switch +fig.fig_state=1; fig.fig_map=0; fig.fig_3d=0; fig.fig_interp=0; disp('Subset simulation domain and convert time') -w=load('w');w=w.w; red=subset_domain(w); disp('Loading and subsetting detections') p=sort_rsac_files(prefix); time_bounds=subset_detection_time(red,p); +red.time_bounds=time_bounds; g = load_subset_detections(prefix,p,red,time_bounds,fig); % Parameters of the objective function @@ -28,10 +30,17 @@ g = load_subset_detections(prefix,p,red,time_bounds,fig); params.alpha=0; params.TC = 1/24; % detection time constants in hours params.stretch=input_num('Tmin,Tmax,Tneg,Tpos',[0.5,10,5,10]); -params.weight=input_num('water,land,low,nominal,high confidence fire',[-1,-1,0.2,0.6,1]); +params.weight=input_num('water,land,low,nominal,high confidence fire',... + [-10,-10,0.2,0.6,1]); params.power=input_num('correction smoothness',1.02); params.doplot=0; params.dx=444; params.dy=444; -J=detection_objective(red.tign,0,g,params); +% objective function is penalty minus log likelihood -> min +J=-detection_objective(red.tign,0,g,params,red); + +s=sprintf('Data log likelihood %18.11e',J); +plot_state(fig.fig_state,red,s,red.tign,g,time_bounds(1:2)); + +end diff --git a/other/Matlab/detection/detection_objective.m b/other/Matlab/detection/detection_objective.m index 023fd3f3..938bb0ca 100644 --- a/other/Matlab/detection/detection_objective.m +++ b/other/Matlab/detection/detection_objective.m @@ -1,7 +1,6 @@ - function varargout=detection_objective(tign,h,g,params) - % [J,delta]=detection_objective(tign,h,g,params) - % J=detection_objective(tign,h,g,params) + function varargout=detection_objective(tign,h,g,params,red) + % [J,delta,f0,f1]=detection_objective(tign,h,g,params) % compute objective function and optionally gradient delta direction % % input: @@ -15,9 +14,14 @@ % to compute data likelihood only, set h=0 and params.alpha=0 % % output - % J objective function - % delta search direction, grad J(tign+h) + % J objective function = penalty minus log likelihood + % delta (optional) preconditioned search direction, grad J(tign+h) + % f0 (optional) contributions of mesh cells to log likelihood + % f1 (optional) contributions of mesh cells to the derivative + % of the log likelihood + like_plots=0; + T=tign+h; f0=0; f1=0; @@ -32,13 +36,18 @@ detections=sum(psi(:)>0); f0=f0+f0k; f1=f1+f1k; - % figure(14);mesh(red.fxlong,red.fxlat,psi),title('psi') - % figure(15);mesh(red.fxlong,red.fxlat,f0k),title('likelihood') - % figure(16);mesh(red.fxlong,red.fxlat,f1k),title('gradient') - % drawnow + if like_plots, + plot_loglike(4,f0k,'Data likelihood',red) + plot_loglike(5,f1k,'Data likelihood derivative',red) + plot_loglike(6,psi,'psi',red) + plot_loglike(7,g(k).time-T,'Time since fire arrival',red) + hold on;contour3(red.fxlong,red.fxlat,g(k).time-T,[0 0],'r') + hold off + drawnow + end end - %figure(15);mesh(red.fxlong,red.fxlat,f0k),title('likelihood') - %figure(16);mesh(red.fxlong,red.fxlat,f1k),title('gradient') + plot_loglike(2,f0,'Data likelihood',red) + plot_loglike(3,f1,'Data likelihood derivative',red) drawnow F=f1; % objective function and preconditioned gradient @@ -48,8 +57,8 @@ J2 = -ssum(f0); J = params.alpha*J1 + J2; fprintf('Objective function J=%g (J1=%g, J2=%g)\n',J,J1,J2); + varargout{1}=J; if nargout==1, - varargout={J}; return end gradJ = params.alpha*Ah + F; @@ -59,9 +68,17 @@ plotstate(8,F,'Detection likelihood derivative',0); plotstate(10,gradJ,'gradient of J',0); end - delta = solve_saddle(params.Constr_ign,h,F,... - @(u) poisson_fft2(u,[params.dx,params.dy],-params.power)/params.alpha); - varargout=[{J},{delta}]; + if params.alpha>0, + delta = solve_saddle(params.Constr_ign,h,F,... + @(u) poisson_fft2(u,[params.dx,params.dy],-params.power)/params.alpha); + varargout{2}=delta; + end + if nargout >= 3, + varargout{3}=f0; + end + if nargout >= 4, + varargout{4}=f1; + end % figure(17);mesh(red.fxlong,red.fxlat,delta),title('delta') % plotstate(11,delta,'Preconditioned gradient',0); %fprintf('norm(grad(J))=%g norm(delta)=%g\n',norm(gradJ,'fro'),norm(delta,'fro')) diff --git a/other/Matlab/detection/plot_loglike.m b/other/Matlab/detection/plot_loglike.m new file mode 100644 index 00000000..55a256ca --- /dev/null +++ b/other/Matlab/detection/plot_loglike.m @@ -0,0 +1,22 @@ +function plot_loglike(fig,f,s,red) +% plot_loglike(fig,f,s,red) +% +% display array in a way suitable for log likelihood +% overay with fireline at end of observations = time_bounds(2) +% fig figure number +% f array to show +% s title +% red the reduced structure with everyhing +figure(fig) +hold off, clf +% show field f +pcolor(red.fxlong,red.fxlat,f) +shading interp +colormap default +colorbar +hold on +% add fireline +contour3(red.fxlong,red.fxlat,red.tign-red.time_bounds(2),[0 0],'k') +hold off +title(s) +end -- 2.11.4.GIT