line_search_onm2

 function  [alpha, img, dv, opt] = line_search_onm2(imgk, dx, data1, img1, N, W, hps2RtR, hpt2LLt, dv0, opt)
 line search function with a fitted polynomial of O(n-2) where n is the number of perturbations
 (C) 2013 Alistair Boyle
 License: GPL version 2 or version 3

0001 function  [alpha, img, dv, opt] = line_search_onm2(imgk, dx, data1, img1, N, W, hps2RtR, hpt2LLt, dv0, opt, retry, pf_max)
0002 % function  [alpha, img, dv, opt] = line_search_onm2(imgk, dx, data1, img1, N, W, hps2RtR, hpt2LLt, dv0, opt)
0003 % line search function with a fitted polynomial of O(n-2) where n is the number of perturbations
0004 % (C) 2013 Alistair Boyle
0005 % License: GPL version 2 or version 3
0006 
0007 perturb = calc_perturb(imgk, dx, opt);
0008 
0009 if nargin < 11
0010   retry = 0;
0011 end
0012 
0013 if nargin < 12
0014   pf_max = length(perturb)-2;
0015 end
0016 
0017 % fwd_solve is the most expensive part generally, count how many we do
0018 if ~isfield(opt, 'fwd_solutions')
0019    opt.fwd_solutions = 0;
0020 end
0021 x = imgk.elem_data;
0022 
0023 if(perturb(1) ~= 0)
0024   error('first perturbation min(inv_model.inv_solve_abs_{GN,CG,core}.line_search.perturb) expects alpha=0');
0025 end
0026 
0027 % Compute the forward model for each perturbation step
0028 img = imgk;
0029 % mlist is our search result for each alpha value, perturb(i)
0030 %  -- NaN: initiailized but not calculated
0031 %  -- -Inf: should not occur we have code that converts calculated NaNs and -Inf to +Inf
0032 %  -- +Inf: calculated value was bad, ignore it
0033 if opt.verbose > 1
0034    fprintf('      i     = ');
0035    fprintf('    [%d]  \t', 1:length(perturb));
0036    fprintf('\n');
0037    fprintf('      alpha = ');
0038    fprintf(' %8.3g\t', perturb);
0039    fprintf('\n');
0040    fprintf('              ');
0041 end
0042 mlist= ones(size(perturb))*NaN; % init
0043 for i = 1:length(perturb)
0044    if (i == 1) && (~isempty(dv0))
0045       % don't bother simulating when alpha=0 (we already have the measurements)
0046       dv = dv0; % @ alpha=0 from the previous line search iteration
0047    else
0048       % fwd_solve and then calculate measurement error (dv)
0049       img.elem_data = x + perturb(i)*dx;
0050       [dv, opt] = feval(opt.line_search_dv_func, img, data1, N, opt);
0051    end
0052    % build de, the change in elem_data from the initial guess
0053    de = feval(opt.line_search_de_func, img, img1, opt);
0054    % we only calculate a new residual if the input data is "sane"
0055    if any(isnan(dv) | isinf(dv))
0056       warning(sprintf('%d of %d elements in dv are NaN or Inf', ...
0057                       length(dv), ...
0058                       length(find(isnan(dv) | isinf(dv)))));
0059       mlist(i) = +Inf;
0060    elseif any(isnan(de) | isinf(de))
0061       warning(sprintf('%d of %d elements in de are NaN or Inf', ...
0062                       length(de), ...
0063                       length(find(isnan(de) | isinf(de)))));
0064       mlist(i) = +Inf;
0065    else
0066       % calculate the residual
0067       mlist(i) = feval(opt.residual_func, dv, de, W, hps2RtR, hpt2LLt);
0068       if any(isnan(mlist(i)) | isinf(mlist(i)))
0069          mlist(i) = +Inf; % NaN or Inf are converted to Inf, since we use NaN to indicate initialized but not calculated
0070       end
0071    end
0072    if opt.verbose > 1
0073       fprintf(' %8.3g\t',mlist(i));
0074    end
0075    if mlist(i)/mlist(1) > 1e10
0076       if opt.verbose > 1
0077          for j=(i+1):length(perturb)
0078             fprintf('   [skip]\t');
0079          end
0080       end
0081       break
0082    end
0083 end
0084 if opt.verbose > 1
0085    fprintf('\n');
0086 end
0087 % drop bad values
0088 if isinf(mlist) % NaN's from any calculation were converted to Inf's earlier
0089    warning('encoutered NaN or +-Inf residuals, something has gone wrong in the line search, converting to large numbers and carrying on');
0090 end
0091 
0092 % TODO looks like this was quiting when there are still good choices left
0093 %if all((mlist/mlist(1)-1) < 1e-4) % < 0.01% change
0094 %   % TODO maybe we need to search *larger* perturbations here... for now we just short circuit the repeated retries at the end, when we are not improving
0095 %   if opt.verbose > 1
0096 %      fprintf('      stopping line search: no further improvements observed\n');
0097 %   end
0098 %   img = imgk;
0099 %   alpha = 0;
0100 %   dv = dv0;
0101 %   return;
0102 %end
0103 
0104 % For our poly fit, we drop all NaN and Inf values
0105 goodi = find((~isnan(mlist)) & (~isinf(mlist)));
0106 alpha=perturb(end);
0107 meas_err = +Inf; % make sure we grab the min(mlist) if we're not doing a polyfit
0108 if length(goodi) > 2
0109   % Select the best fitting step, we scale and
0110   p_rng = range(perturb(goodi)); % p_min = 0
0111   pfx = log10(perturb(goodi)/p_rng);
0112   pfx(1) = -100; % log10(0) = -Inf --> -1e100 so that it's finite
0113   pf= polyfit(pfx, mlist(goodi), length(goodi)-2);
0114   % search for the function minima in the range perturb(2:end)
0115   %   pf(1)*log10(x).^2+pf(2)*log10(x)+pf(3);
0116   FF = @(pf, x) polyval(pf, log10(x/p_rng));
0117   alpha = fminbnd(@(x) FF(pf, x), perturb(min(goodi)), perturb(end));
0118   % now check how we did
0119   img.elem_data = x + alpha*dx;
0120   [dv, opt] = feval(opt.line_search_dv_func, img, data1, N, opt);
0121   de = feval(opt.line_search_de_func, img, img1, opt);
0122   meas_err = feval(opt.residual_func, dv, de, W, hps2RtR, hpt2LLt);
0123   if opt.verbose > 1
0124      fprintf('      step size = %0.3g, misfit = %0.3g, expected = %0.3g\n', alpha, meas_err, FF(pf, alpha));
0125   end
0126 
0127   % check how close we were to the line fit estimate
0128   % suggest stronger regularization if we're way off
0129   % (we could automate this update if we wanted to turn on some hueristic behaviour)
0130   est_err = meas_err / FF(pf, alpha);
0131   if (opt.verbose > 1) && ((est_err > 1.3) || (est_err < 0.5))
0132     fprintf('      step misfit missed estimate (%0.1fx est)\n', est_err);
0133     fprintf('        consider stronger regularization?\n');
0134   end
0135 else % only two points
0136   % we provide a FF and pf that will work for plot_line_optimize()
0137   % this is a straight line between alpha=1 and alpha=1/10
0138   pf = [];
0139   FF = @(pf, x)  -(mlist(1) - mlist(end))*0.8*log10(x) + mlist(end);
0140 end
0141 % We save our first choice, in case we are plotting the line search
0142 alpha1 = alpha; % better guess: minima of the fitted curve
0143 meas_err1 = meas_err;
0144 
0145 % what if we're making a bad choice?
0146 % if our choice sucked, we've already calculated mlist(:) other choices, go with the minimum
0147 if meas_err > min(mlist)
0148   [meas_err,mi]= min(mlist);
0149   alpha= perturb(mi);
0150   img.elem_data = x + alpha*dx;
0151   dv = []; % we'll need to recalculate this later since we didn't keep it
0152   if (length(goodi) > 2) && (opt.verbose > 1)
0153     fprintf('      did not like our step selection - choose one of perturb values\n');
0154   end
0155 end
0156 
0157 if opt.verbose > 1
0158    max_alpha_str = '';
0159    if alpha > perturb(end)-eps
0160      max_alpha_str = ' (max)';
0161    end
0162    fprintf('      step size = %0.3g%s, misfit = %0.3g selected\n', alpha, max_alpha_str, meas_err);
0163 end
0164 %
0165 % keyboard
0166 % must create plots before changing the perturb values
0167 if opt.line_search_args.plot
0168    clf;
0169    plot_line_optimize(perturb, mlist, alpha, meas_err, alpha1, meas_err1, FF, pf);
0170    if isfield(opt,'fig_prefix') % TODO assign from base options if set
0171       k=1; % TODO iteration count; TODO retry count
0172       print('-dpdf',sprintf('%s-ls%d-retry%d',opt.fig_prefix,k,retry));
0173       print('-dpng',sprintf('%s-ls%d-retry%d',opt.fig_prefix,k,retry));
0174       saveas(gcf,sprintf('%s-ls%d-retry%d.fig',opt.fig_prefix,k,retry));
0175    end
0176 end
0177 
0178 % update perturbations
0179 if meas_err >= mlist(1)
0180     if mlist(1)*1.05 < mlist(goodi(end))
0181        % this happens when the solution blew up -- the measurement fit was worse than if we did nothing
0182        if opt.verbose > 1
0183           fprintf('      reducing perturbations /10: bad step\n');
0184        end
0185        % try a smaller step next time (10x smaller)
0186        % this keeps the log-space distance between sample points
0187        perturb = perturb/10;
0188     elseif perturb(end) > 1.0-10*eps
0189        if opt.verbose > 1
0190           fprintf('      expanding perturbations x10: ... but we''d be searching past alpha=1.0, giving up\n');
0191        end
0192        return % we give up early
0193     elseif perturb(end)*10 > 1.0-10*eps
0194        if opt.verbose > 1
0195           fprintf('      expanding perturbations (limit alpha=1.0): bad step\n');
0196        end
0197        perturb = perturb/perturb(end); % ... max(alpha)=1.0
0198     else % we didn't really get any difference in solutions
0199        % this happens when the perturbations are too small, we are too close to
0200        % the current solution
0201        if opt.verbose > 1
0202           fprintf('      expanding perturbations x10: bad step\n');
0203        end
0204        % try a larger step next time (10x larger)
0205        % this keeps the log-space distance between sample points
0206        perturb = perturb*10;
0207     end
0208 else % good step
0209     % stretch out the perturbations if we're not making much progress
0210     if all(mlist(goodi)/mlist(1)-1 > -10*opt.dtol) && ...
0211        (perturb(end)*10 < 1.0+10*eps)
0212        if opt.verbose > 1
0213           fprintf('      expand perturbations x10 for next iteration\n');
0214           fprintf('      (didn''t make much progress this iteration)\n');
0215        end
0216        opt.line_search_args.perturb = opt.line_search_args.perturb*10;
0217     else % or just recentre around our best answer
0218        % this keeps the log-space distance between sample points but
0219        % re-centres around the most recent alpha
0220        if opt.verbose > 1
0221           fprintf('      update perturbations around step = %0.3g (limit alpha=1.0)\n', alpha);
0222        end
0223        if alpha/perturb(end)*2 > 1.0 - 10*eps
0224           perturb = perturb/perturb(end);
0225        else
0226           perturb = perturb*(alpha/perturb(end))*2;
0227        end
0228     end
0229 end
0230 % jiggle the perturb values by 1% --> if we're stuck in a recursion
0231 % of bad perturb values maybe this is enough to break us out
0232 perturb = perturb .* exp(randn(size(perturb))*0.01);
0233 % fix if we exceeded alpha=1.0
0234 if perturb(end) > 1.0 - eps
0235    perturb = perturb/perturb(end);
0236 end
0237 opt.line_search_args.perturb = perturb;
0238 
0239 % Record the corresponding parameters
0240 %img.elem_data= exp(img.logCond);
0241 
0242 % we took a bad step, try again but don't recurse indefinitely
0243 if alpha == 0 && retry < 5
0244   if opt.verbose > 1
0245      fprintf('    retry#%d (attempt with new perturbations)\n', retry+1);
0246   end
0247   [alpha, img, dv, opt] = line_search_onm2(imgk, dx, data1, img1, N, W, hps2RtR, hpt2LLt, dv0, opt, retry+1, pf_max);
0248 end
0249 
0250 %%%
0251 % calculate search values for \alpha
0252 % 1. sort in ascending order
0253 % 2. scale to a range within finite precision (don't waste time searching when
0254 %    results will be inf/nan)
0255 function perturb = calc_perturb(imgk, dx_in, opt)
0256 if opt.verbose > 1
0257    disp('      line search (finite precision) limits');
0258 end
0259 % scale
0260 % When log numbers are converted to base numbers, they frequently result in Inf
0261 % when dx is large.  We scale 'perturb' here, so that we are line searching
0262 % within a numerically stable region for our finite precision numbers.
0263 %   log(realmax) = 709.7827
0264 % log10(realmax) = 308.2547
0265 % - canonicalize the img data, so we don't have to deal with default forms
0266 if ~isfield(imgk, 'current_params')
0267    imgk.current_params = {'conductivity'};
0268 end
0269 if ~isfield(imgk, 'params_sel')
0270    imgk.params_sel = {[1:length(imgk.elem_data)]};
0271 end
0272 if ~iscell(imgk.current_params)
0273    imgk.current_params = {imgk.current_params};
0274 end
0275 if ~iscell(imgk.params_sel)
0276    imgk.params_sel = {imgk.params_sel};
0277 end
0278 
0279 if isfield(imgk, 'inv_model') && isfield(imgk.inv_model, 'fwd_model')
0280    md = max(range(imgk.inv_model.fwd_model.nodes)); % model range (coordinates)
0281 end
0282 
0283 % determine the maximum alpha
0284 max_alpha = +inf;
0285 min_alpha = +inf;
0286 for i=1:length(imgk.current_params)
0287    p = imgk.current_params{i};
0288    s = imgk.params_sel{i};
0289    x = imgk.elem_data(s);
0290    dx = dx_in(s);
0291    % TODO these could be based on the limits provided as args to inv_solve_abs_GN, instead they are hardcoded here...
0292    is_mvmt = (length(p) >= 8) && strcmp(p(end-7:end),'movement');
0293    if strcmp(p(1:4), 'log_')
0294       lp = log(realmax/2); % largest positive floating point number (double): Limit_Positive
0295       ln = -inf; % or = log(realmin/2); % largest negative floating point number (double): Limit Negative
0296       % for log space, we should have an ln = -inf --> exp(-900) = 0
0297       if is_mvmt
0298          lp = log(md);
0299       end
0300    elseif strcmp(p(1:6), 'log10_')
0301       lp = log10(realmax/2); % largest positive floating point number (double): Limit_Positive
0302       ln = -inf; % or = log10(realmin/2); % largest negative floating point number (double): Limit Negative
0303       % for log10 space, we should have an ln = -inf --> 10.^-900 = 0
0304       if is_mvmt
0305          lp = log10(md);
0306       end
0307    else
0308       lp = +realmax/2;
0309       ln = -realmax/2;
0310       if is_mvmt
0311          lp = +md;
0312          lp = -md;
0313       end
0314    end
0315    % lower limit on \alpha prior to x = x + alpha*dx --> +/-inf; % (explodes)
0316    %   \alpha_min = ((max or min) - x) / \delta_x
0317    au=(lp-x)./dx; au(dx<=0)=NaN; au(isnan(au))=+inf; au=min(au);
0318    a_max = au;
0319    au=(ln-x)./dx; au(dx>=0)=NaN; au(isnan(au))=+inf; au=min(au);
0320    if (au < a_max)
0321       a_max = au;
0322    end
0323    % lower limit on \alpha prior to x == x + alpha*dx; % (no change)
0324    %   \alpha_min = \epsilon / \delta_x
0325    if is_mvmt
0326       al=1e-3./abs(dx); % don't care about movement less than 1mm
0327       % TODO configurable? 'reconstruction tolerance'?
0328    else
0329       al=eps(x)./abs(dx);
0330    end
0331    al(isinf(al))=NaN; al(isnan(al))=realmax; al=min(al);
0332    if isnan(al)
0333       a_min = 0;
0334    else
0335       a_min = al;
0336    end
0337    if opt.verbose > 1
0338       fprintf('        %s: alpha range = %0.3g -- %0.3g\n', p, a_min, a_max);
0339    end
0340    % adjust global limits
0341    if a_max < max_alpha
0342       max_alpha = a_max;
0343    end
0344    if a_min < min_alpha
0345       min_alpha = a_min;
0346    end
0347 end
0348 
0349 % sort
0350 p=sort(opt.line_search_args.perturb);
0351 % scale
0352 if (p(end) > max_alpha) || (p(2) < min_alpha)
0353    p(p<realmin/2) = [];
0354    p=log10(p); ap=log10(max_alpha); an=log10(min_alpha);
0355    if range(p) >  ap-an
0356       p=p*(ap-an)/range(p);
0357    end
0358    if p(end) > ap
0359       p=p-(max(p)-ap);
0360    elseif p(1) < an
0361       p=p+(an-min(p));
0362    end
0363    p=[0 10.^p];
0364    if opt.verbose > 1
0365       fprintf('        alpha (before) = ');
0366       fprintf('%0.3g ', sort(opt.line_search_args.perturb));
0367       fprintf('\n');
0368       fprintf('        alpha (after)  = ');
0369       fprintf('%0.3g ', p);
0370       fprintf('\n');
0371    end
0372 else
0373    if opt.verbose > 1
0374       fprintf('        alpha (unchanged) = ');
0375       fprintf('%0.3g ', p);
0376       fprintf('\n');
0377    end
0378 end
0379 perturb=p;
0380 
0381 
0382 %%% plot the line optimization results
0383 % 1. search locations
0384 % 2. line fit
0385 % 3. selected minima and test point result
0386 % 4. selected \alpha
0387 function plot_line_optimize(perturb, mlist, alpha, meas_err, alpha1, meas_err1, FF, pf)
0388 semilogx(perturb(2:end),mlist(2:end),'xk', 'MarkerSize',10);
0389 hold on;
0390 semilogx(alpha, meas_err,'or', 'MarkerSize',10);
0391 semilogx(alpha1, FF(pf, alpha1), 'ob', 'MarkerSize',10);
0392 semilogx(alpha1, meas_err1, 'xb', 'MarkerSize',10);
0393 % construct the fitted line for plotting
0394 if isnan(perturb(2))
0395    perturb(2) = perturb(end)/10;
0396 end
0397 p= logspace(log10(perturb(2)),log10(perturb(end)),50);
0398 semilogx(p, FF(pf, p),'k','linewidth',2);
0399 semilogx(p, p*0+mlist(1),'k--','linewidth',1); % alpha=0 should be plotted as well
0400 %% *&^%@!# Idiot matlab now doesn't allow "selected" in legend!
0401 legend('perturb', 'selected ', '1st est', '1st act', 'fit', '\alpha=0');
0402 legend('Location', 'EastOutside');
0403 m = [mlist meas_err meas_err1];
0404 mi=find(isnan(m) | isinf(m)); m(mi) = []; % remove bad values
0405 mr = range(m);
0406 if mr < max(m)*1e-14
0407    mr = 1e-14;
0408 end
0409 axis([perturb(2) perturb(end) min(m)-mr*0.2 max(m)+mr*0.2]);
0410 xlabel('step size \alpha'); %,'fontsize',20,'fontname','Times')
0411 ylabel('normalized residuals'); %,'fontsize',20,'fontname','Times')
0412 title({sprintf('best alpha = %1.2e',alpha), ...
0413        sprintf('norm w/o step = %0.4e',mlist(1))}); %,'fontsize',30,'fontname','Times')
0414 %set(gca,'fontsize',20,'fontname','Times');
0415 drawnow;
0416 
0417 function x=range(y)
0418 x=max(y)-min(y);