detection/detection_objective.m

   1     function varargout=detection_objective(tign,h,g,params,red,varargin)
   2     % J              =detection_objective(tign,h,g,params,red)
   3     % [J,delta]      =detection_objective(tign,h,g,params,red)
   4     % [J,delta,f0,f1]=detection_objective(tign,h,g,params,red)
   5     % compute objective function and optionally gradient delta direction
   6     %
   7     % input:
   8     %     tign   prior state
   9     %     h      increment
  10     %     g      detection structur?e from load_subset_detections
  11     %     params parameters set in detect_fit_level2
  12     %
  13     % The objective function is
  14     % J = h'*A*h - log likelihood of tign+h
  15     % to compute data likelihood only, set h=0 and params.alpha=0
  16     %
  17     % output
  18     %     J      objective function = penalty minus log likelihood
  19     %     delta  (optional) preconditioned search direction, grad J(tign+h)
  20     %     f0     (optional) contributions of mesh cells to log likelihood
  21     %     f1     (optional) contributions of mesh cells to the derivative
  22     %                   of the log likelihood
  23
  24     if nargin >=6
  25         new_like = 1;
  26         p_like_spline = varargin{1};
  27         p_deriv_spline = varargin{2};
  28         n_deriv_spline = varargin{3};
  29     else
  30         new_like = 0;
  31     end
  32
  33     like_plots=0;
  34
  35         T=tign+h;
  36         f0=0;
  37         f1=0;
  38         %%% altered weighting of the different data sources, store params
  39         %%% first
  40         params_bak = params;
  41
  42         %% different weighting scheme if perimter data present
  43         perims_present = 0;
  44         for k = 1:length(g)
  45             if strcmp(g(k).file(1:3),'PER')
  46                 if ~perims_present
  47                     fprintf('Perimeter data present \n');
  48                     perims_present = input_num('1 to adjust satellite weights ', 1);
  49                 end
  50             end
  51         end
  52
  53         % compute new likelihood stuff here,
  54
  55         if new_like
  56             fprintf('Using likelihood splines.\n')
  57             % load splines.mat;
  58             %[like_spline,fitness1] = createFit(t,log_like);
  59             %[deriv_spline,fitness2] = createFit(t,log_like_deriv);
  60         end
  61
  62         for k=1:length(g)
  63             % load params in case it has been changed
  64             params = params_bak;
  65
  66 % routine for weighting modis and virrs less when [perimeter data present
  67             non_perim_weight = 0.5;
  68             if ~strcmp(g(k).file(1:3),'PER') && perims_present
  69                 fprintf('Changing weight of satellite detection \n');
  70                 params.weight(3:5) = non_perim_weight*params.weight(3:5);
  71             end
  72
  73             %weight by time
  74             time_weight = 0;
  75             if time_weight
  76                 pwr = 0.5;
  77                 %multiplier = exp((g(k).time-g(length(g)).time))^pwr;
  78                 multiplier = ((g(k).time-min(tign(:)))/(g(length(g)).time-min(tign(:))))^pwr;
  79                 params.weight = multiplier*params.weight;
  80                 fprintf('Changing weights by %f \n',multiplier);
  81             end
  82
  83 % routine for changing modis weights
  84 %             modis_weight = 0.2;
  85 %             if strcmp(g(k).file(1:3),'MOD') | strcmp(g(k).file(1:3),'MYD')
  86 %             %if g(k).file(1:3) == 'MOD' | g(k).file(1:3) == 'MYD'
  87 %                 fprintf('Changing detection weights for MODIS data \n');
  88 %                 params.weight(3:5) = modis_weight*params.weight(3:5);
  89 %             end
  90             psi = ...
  91                 + params.weight(1)*(g(k).fxdata==3)...
  92                 + params.weight(2)*(g(k).fxdata==5)...
  93                 + params.weight(3)*(g(k).fxdata==7)...
  94                 + params.weight(4)*(g(k).fxdata==8)...
  95                 + params.weight(5)*(g(k).fxdata==9);
  96             if new_like
  97                 [f0k,f1k] = evaluate_likes(psi,g(k).time-T,p_like_spline,p_deriv_spline,n_deriv_spline);
  98                 f0k = f0k - max(f0k(:));
  99                 f1k = f1k-max(f1k(:));
 100                 %f1k = -f1k; %% why
 101             else
 102                 [f0k,f1k]=like2(psi,g(k).time-T,params.TC*params.stretch);
 103             end
 104             detections=sum(psi(:)>0);
 105             f0=f0+f0k;
 106             f1=f1+f1k;
 107 %             if k == length(g)
 108 %                 figure,mesh(f0),title('likelihood');
 109 %                 figure,mesh(f1),title('derivative');
 110             if like_plots>1,
 111                 plot_loglike(4,f0k,'Data likelihood',red)
 112                 plot_loglike(5,f1k,'Data likelihood derivative',red)
 113                 plot_loglike(6,psi,'psi',red)
 114                 plot_loglike(7,g(k).time-T,'Time since fire arrival',red)
 115                 hold on;contour3(red.fxlong,red.fxlat,g(k).time-T,[0 0],'r')
 116                 hold off
 117                 drawnow
 118             end
 119         end
 120         if like_plots>0,
 121             plot_loglike(2,f0,'Data likelihood',red)
 122             plot_loglike(3,f1,'Data likelihood derivative',red)
 123         drawnow
 124         end
 125         F=f1;
 126         % objective function and preconditioned gradient
 127         Ah = poisson_fft2(h,[params.dx,params.dy],params.power);
 128         % compute both parts of the objective function and compare
 129         J1 = 0.5*(h(:)'*Ah(:));
 130         J2 = -ssum(f0);
 131         J = params.alpha*J1 + J2;
 132         fprintf('Objective function J=%g (J1=%g, J2=%g)\n',J,J1,J2);
 133         varargout{1}=J;
 134         if nargout==1,
 135             return
 136         end
 137         gradJ = params.alpha*Ah + F;
 138         fprintf('Gradient: norm Ah %g norm F %g\n', norm(Ah,2), norm(F,2));
 139         if params.doplot,
 140             plotstate(7,f0,'Detection likelihood',0);
 141             plotstate(8,F,'Detection likelihood derivative',0);
 142             plotstate(10,gradJ,'gradient of J',0);
 143         end
 144         if params.alpha>0,
 145             delta = solve_saddle(params.Constr_ign,h,F,0,...
 146                 @(u) poisson_fft2(u,[params.dx,params.dy],-params.power)/params.alpha);
 147             bump_remove = 1;
 148             if bump_remove
 149                 % change search to keep bump from forming here
 150                 temp_h = -delta/big(delta);
 151                 temp_analysis = tign+temp_h;
 152                 %figure,mesh(temp_analysis)
 153                 corner_time = min([temp_analysis(1,1),temp_analysis(1,end),...
 154                     temp_analysis(end,1),temp_analysis(end,end)]);
 155                 if max(temp_analysis(:)) > corner_time
 156                     fprintf('Bump forming, inside the detection_objective function \n');
 157                     %figure,mesh(temp_analysis);
 158                     mask_h = temp_analysis >= corner_time;
 159                     delta(mask_h) = 0;
 160                     %figure,mesh(delta);
 161                 end
 162             end %bump removal
 163             varargout{2}=delta; %search in detect_fit_level2
 164         end
 165         if nargout >= 3,
 166             varargout{3}=f0;
 167         end
 168         if nargout >= 4,
 169             varargout{4}=f1;
 170         end
 171         % figure(17);mesh(red.fxlong,red.fxlat,delta),title('delta')
 172         % plotstate(11,delta,'Preconditioned gradient',0);
 173         %fprintf('norm(grad(J))=%g norm(delta)=%g\n',norm(gradJ,'fro'),norm(delta,'fro'))
 174     end
 175
 176