% This codes solves the Generalized Trust Region Subproblem:
%
% min x'Ax - 2a'x   s.t.   lb <= x'Bx <= ub
% where B is indefinite nonsingular, A is positive definite
%
% The code starts with preprocessing that either solves the problem by 
% returning an interior solution, an explicit solution in hard case 2
% instances via shift and deflation, or reduces the problem to an
% equivalent equality constrained problem, i.e., x'Bx = s.
%
% The equality constrained problem is further reformulated into
%                  max k(t) := k0(t) - t
% where k0(t) = inf{y'[t -a';-a A]y: y'[1 0;0 B]y = s + 1},
% and solved via the extended Rendl Wolkowicz algorithm.
%
% Input
% 
% A, B, a, lb, ub - input data from the problem
% gtol            - termination tolerance
% t0, lambda0     - t0 is inf{lambda + a'(A - lambda B)^{-1}a: underline(lambda) < lambda < overline(lambda)}, 
%                   where underline(lambda) and overline(lambda)} are the largest negative and smallest positive 
%                   generalized eigenvalue of the pencil (A,B), respectively. 
%                   lambda0 is the corresponding value of lambda corresponding to this t0.
%
% Output (argout in order)
%
% y       - approximate optimal solution
% iter    - number of iterations to solve the equality constrained problem;
% fval    - approximate primal optimal value: y'Ay - 2a'y
% feas    - approximate primal feasibility: y'By - min(max(y'By,lb),ub)
% dval    - approximate dual optimal value: k(t)
% k_slope - slope of k at termination
% t_final - total CPU time
% lambda  - approximate optimal Lagrange multiplier. Equals NaN if
%           initialization fails.
% s       - the equality constraint that is attained. Equals NaN for interior
%           solution.
%
% Subroutines called: Ahandle, Bhandle, Chandle, Dhandle, eigifp
%
% Written by Ting Kei Pong and Henry Wolkowicz
% Last updated on Nov 9, 2012

function varargout = gtrsalgo_ul_indB(A,a,B,lb,ub,gtol,t0,lambda0)

tstart = tic;
global A0
global B0
global a0
global tt
global s
global shift
global flag_shift_lambda
global MVcount

A0 = A;
B0 = B;
a0 = a;
MVcount = 0;

fprintf(' ************** START ***************\n')

nA = normest(A,1);
na = norm(a);

fprintf(' ------------- Initialization -------------\n')
zerotol = 1e-8;
[s,flag_case,output] = GTRS_case_check(A,a,B,lb,ub,zerotol);

if flag_case ~= 0 %  interior solution or hardcase 2
  t_final = toc(tstart);
  
  fprintf('\n ------------- Results -------------\n')

  y = output.y;
  fval = output.fval;
  feas = output.feas;
  lopt = output.lopt;
  fprintf(' Primal value       %16.10f  Feasibility      %16.10f\n', fval, feas)
  fprintf(' Rel Opt Cond.      %16.10f  CPU time         %16.10f\n', norm(A*y - a - lopt*(B*y))/(nA + na + 1), t_final)
  fprintf(' Number of MV       %16.10f\n', MVcount)
  
  if nargout >= 1; varargout{1} = y; end
  if nargout >= 2; varargout{2} = 0; end
  if nargout >= 3; varargout{3} = fval; end
  if nargout >= 4; varargout{4} = feas; end
  if nargout >= 5; varargout{5} = fval; end
  if nargout >= 6; varargout{6} = 0; end
  if nargout >= 7; varargout{7} = t_final; end
  if nargout >= 8; varargout{8} = lopt; end
  if nargout >= 9; varargout{9} = s; end
  return
end

n = size(A,1);
s1 = abs(s + 1);
ss1 = sqrt(s1);
fprintf(' Problem size: %d, s: %16.10f\n', n,s);

tol = 1e-14; % for pcg
eigtol = zerotol;

if nargin <= 5
  % check assumption: positive slope at t(0)
  [u,dummy] = pcg(Ahandle,a,tol,1000);
  stmp = u'*B*u;
  if s <= stmp || stmp <= -1 % sufficient condition only
    fprintf(' Initialization fails\n')
    if nargout >= 1; varargout{1} = 0; end
    if nargout >= 2; varargout{2} = 0; end
    if nargout >= 3; varargout{3} = 0; end
    if nargout >= 4; varargout{4} = 0; end
    if nargout >= 5; varargout{5} = 0; end
    if nargout >= 6; varargout{6} = 0; end
    if nargout >= 7; varargout{7} = 0; end
    if nargout >= 8; varargout{8} = 0/0; end
    if nargout >= 9; varargout{9} = s; end
    return
  end
  low = a'*u;
  shift = 0;
  flag_shift_lambda = 0;
else
  low = t0 + 1; % heuristic to get left end point
  shift = lambda0; % so that Lt is positive definite...
  flag_shift_lambda = 1;
  divide = 0;
end

opts.initialvec = randn(n+1,1);
opts.tol = eigtol;
opts.disp = 0;
opts.useprecon = 'NO';
opts.maxit = 10000;
opts.size = n+1;

tinit2 = tic;
low = low - shift; % shift of A is performed in Chandle
tt = low;
[v, bv_init, dummy] = GTRS_ktsubproblem;
fprintf(' eig diff %g lambda %g divide %g\n', norm(Dhandle(v) - dummy*Chandle(v)), dummy, divide)

neg_slope = 0; % record existence of point with negative slope

k_slope = s1*v(1)^2 - 1;
% backtrack to get left end point
while k_slope < 0 && divide < 4;  
  high = tt;
  dval = s1*(-1/dummy) - tt;
  ptgood = [high abs(v(1)) k_slope dval];
  gv = v;  
  neg_slope = 1;
  
  divide = divide + 1;
  low = t0 + 1/5^divide; % heuristic to get left end point
  low = low - shift; % shift of A is performed in Chandle
  tt = low;
  opts.initialvec = bv_init; % warm start
  [v, bv_init, dummy] = GTRS_ktsubproblem;
  fprintf(' eig diff %g lambda %g divide %g\n', norm(Dhandle(v) - dummy*Chandle(v)), dummy, divide)
  k_slope = s1*v(1)^2 - 1;
end
if divide >= 4
  error(' initialization fails')
end
dval = s1*(-1/dummy) - tt;
ptbad = [low abs(v(1)) k_slope dval]; % t v(1) slope
bv = v;
gv_init = bv_init;

if neg_slope == 0 % no point with negative slope available  
  while 1 == 1 % loop to find right end point
    tt = 1.2*tt; % heuristic; tt must be positive
    opts.initialvec = gv_init;
    [v, gv_init, dummy] = GTRS_ktsubproblem;
    fprintf(' eig diff %g lambda %g\n', norm(Dhandle(v) - dummy*Chandle(v)), dummy)
    k_slope = s1*v(1)^2 - 1;
    fprintf(' finding negative slope\n')
    if k_slope > 0
      low = tt; % update left end point info
      dval = s1*(-1/dummy) - tt;
      ptbad = [low abs(v(1)) k_slope dval]; % t v(1) slope
      bv = v;
    else
      break
    end
  end
  high = tt;
  dval = s1*(-1/dummy) - tt;
  ptgood = [high abs(v(1)) k_slope dval];
  gv = v;
end

fprintf(' Initial end points: %16.10f  %16.10f; time %g\n', low, high, toc(tinit2))

side = 0;
iter = 0;

gap = inf;
fbest = inf;
feas = inf;
optcon = inf;
k_slope = 1;

% Enter main loop
fprintf(' ------------- Enter main loop -------------\n')
v_init = gv_init;

while ~((abs(k_slope)/s1 < sqrt(gtol) && abs(gap)/(abs(fbest) + 1) < gtol && abs(feas)/(abs(s) + 1) < gtol && optcon < sqrt(gtol)) || (high - low)/(abs(high) + abs(low)) < gtol/100) && iter < 30
    
  iter = iter + 1;
  
  gindex = size(ptgood,1);
  bindex = size(ptbad,1);
  
  % strategies for maximizing k
  tt = (high + low)/2; % midpoint rule  
  
  [tt,high,low] = vertical_cut(gindex,bindex,ptgood,ptbad,high,low,tt); % vertical cut
  
  tt = trig_interpol(gindex,bindex,ptgood,ptbad,high,low,tt); % triangle interpolation; replaces the previous only when triangle lies between low and high
  
  tt = inv_interpol(gindex,bindex,ptgood,ptbad,side,ss1,high,low,tt); % inverse interpolation; replaces the previous only when inv lies between low and high

  eigtol = max(min(abs(k_slope)*1e-4,(high - low)/(abs(high) + abs(low))),zerotol); % approximate eig % has to be fairly accurate
%   eigtol = zerotol; % approx eig fails when abs(s) is large
  if side > 0
    v_init = gv_init;
  elseif side < 0
    v_init = bv_init;
  end

  opts.initialvec = v_init;
  opts.tol = eigtol;
  opts.maxit = 5000;
  [v, v_init, dummy] = GTRS_ktsubproblem;
  vCv = 1/(-dummy);
%   fprintf(' eig diff %g lambda %g\n', norm(Dhandle(v) - dummy*Chandle(v)), dummy)
  
  k_slope = s1*v(1)^2 - 1;
  dval = s1*vCv - tt;
  
  % record good side/bad side points
  
  if k_slope > 0
    
    ptbad = [tt abs(v(1)) k_slope dval;ptbad];
    low = tt;
    side = min(-1,side - 1);
    bv = v;
    bv_init = v_init;
    
    if gindex > 0 % primal step to the boundary

      vf = v(2:end)/v(1);
      if abs(gv(1)) < zerotol 
        z = gv(2:end); % take a step with smaller magnitude
        fprintf(' warning: first entry small \n')
      else
        z = gv(2:end)/gv(1) - vf; % take a step between 0 and 1
      end
      Bz = Bhandle(z);
      zBz = z'*Bz;
      xyprod = Bz'*vf;  % inner product
      ndiff = s - vf'*Bhandle(vf);
      sxy = sign(xyprod);
      if sxy == 0,
        sxy = 1;
      end
      tau = ndiff/(xyprod + sxy*sqrt(xyprod^2 + zBz*ndiff));
      vnew = vf + tau*z;
      fval = vnew'*Ahandle(vnew) - 2*a'*vnew - shift*s;
      if fbest >= fval
        fbest = fval;
        y = vnew;
      end
      gap = fbest - dval;
      Bytmp = Bhandle(y);
      feas = y'*Bytmp - s;
      lambdaout = sign(s + 1)*vCv + shift;
      optcon = norm(Ahandle(y) - lambdaout*Bytmp - a)/(nA + na + 1);
      fprintf('%3d bad : shifted primal value %6.5e, gap %3.2e, feas %3.2e\n', iter, fval, gap, feas) 
      
    end
    
  elseif k_slope < 0
    
    ptgood = [tt abs(v(1)) k_slope dval;ptgood];
    high = tt;
    side = max(1,side + 1);
    gv = v;
    gv_init = v_init;
    
    if bindex > 0 % primal step to the boundary
      
      vf = bv(2:end)/bv(1);
      if abs(v(1)) < zerotol
        z = v(2:end); % take a step with smaller magnitude
        fprintf(' warning: first entry small \n')
      else
        z = v(2:end)/v(1) - vf; % take a step between 0 and 1
      end
      Bz = Bhandle(z);
      zBz = z'*Bz;
      xyprod = Bz'*vf;  % inner product
      ndiff = s - vf'*Bhandle(vf);
      sxy = sign(xyprod);
      if sxy == 0,
        sxy = 1;
      end
      tau = ndiff/(xyprod + sxy*sqrt(xyprod^2 + ndiff*zBz));
      vnew = vf + tau*z;
      fval = vnew'*Ahandle(vnew) - 2*a'*vnew - shift*s;
      if fbest >= fval
        fbest = fval;
        y = vnew;
      end
      gap = fbest - dval;
      Bytmp = Bhandle(y);
      feas = y'*Bytmp - s;
      lambdaout = sign(s + 1)*vCv + shift;
      optcon = norm(Ahandle(y) - lambdaout*Bytmp - a)/(nA + na + 1);
      fprintf('%3d good: shifted primal value %6.5e, gap %3.2e, feas %3.2e\n', iter, fval, gap, feas) 
      
    end
    
  end

end

dval = dval + shift*s;
lambda = sign(s + 1)*vCv + shift;
fval = y'*Ahandle(y) - 2*a'*y;
constraint = y'*Bhandle(y);
feas = constraint - max(min(constraint,ub),lb);

t_final = toc(tstart);

% print output

fprintf('\n ------------- Results -------------\n')

fprintf(' Primal value       %16.10f  Dual value       %16.10f\n', fval, dval)
fprintf(' Gap                %16.10f  Feasibility      %16.10f\n', gap, feas)
fprintf(' interval width     %16.10f  Rel Opt Cond.    %16.10f\n', high - low, optcon)
fprintf(' k_slope (approx)   %16.10f  CPU time         %16.10f\n', k_slope, t_final)
fprintf(' Number of MV       %16.10f\n', MVcount)

if nargout >= 1; varargout{1} = y; end
if nargout >= 2; varargout{2} = iter; end
if nargout >= 3; varargout{3} = fval; end
if nargout >= 4; varargout{4} = feas; end
if nargout >= 5; varargout{5} = dval; end
if nargout >= 6; varargout{6} = k_slope; end
if nargout >= 7; varargout{7} = t_final; end
if nargout >= 8; varargout{8} = lambda; end
if nargout >= 9; varargout{9} = s; end





% looped functions

function [s,flag_case,output] = GTRS_case_check(A,a,B,lb,ub,zerotol)

global flag_shiftA
global vall
global alphall 
global lambdalow
global lambdaup

flag_shiftA = 0;          % flag_shiftA will be changed to 1 if we have the possible hard case
n = size(A,1);

opts.useprecon = 'NO';
opts.disp = 0;
opts.maxit = 10000;
opts.initialvec = sum(A)';
opts.tol = zerotol;
opts.size = n;

tinit1 = tic;
% Calculate the first positive and negative eigenvalues, with the corresponding eigenvectors, of the matrix pencil (A,B)
[mu,v3] = eigifp('nBhandle','Ahandle',1,opts);
lambdaup = -1/mu;
[mu,v1] = eigifp('Bhandle','Ahandle',1,opts);
lambdalow = 1/mu;
fprintf(' time for finding interval: %g; interval %g, %g\n', toc(tinit1), lambdalow, lambdaup)

if lambdalow*lambdaup >= 0
  error(' Initialization failed: cannot locate initial interval.\n')
end

valow = v1/norm(v1);
vaup = v3/norm(v3);
tol = 1e-14; % for pcg

vaaup = vaup'*a/na;
vaalow = valow'*a/na;
if vaaup < 0 % loop to ensure the inner product is positive
  vaup = -vaup;
  vaaup = -vaaup;
end
if vaalow < 0 % loop to ensure the inner product is positive
  valow = -valow;
  vaalow = -vaalow;
end
fprintf(' vaaup %g vaalow %g\n', vaaup, vaalow)
    
% solve for possible interior solution
[xtemp, dummy] = pcg('Ahandle',a,tol,1000);  
fprintf(' pcg residual %g\n', norm(Ahandle(xtemp) - a))
Bxtemp = Bhandle(xtemp);
tau1 = xtemp'*Bxtemp - lb;
tau2 = ub - xtemp'*Bxtemp;

if tau1 >=0 && tau2 >= 0 % interior solution found
      
  fprintf('** Interior solution\n');
  xopt = xtemp;
  lopt = 0;    
  fopt = xopt'*Ahandle(xopt) - 2*a'*xopt;
               
  flag_case = 2; % interior solution
  output.y = xopt;
  output.fval = fopt;
  output.lopt = lopt;
  constraint = xopt'*Bhandle(xopt);
  output.feas = constraint - max(min(constraint,ub),lb);
  s = 0/0;
  return            
end
    
if tau1 < 0; flag_kcase = 0; end % maximize k(t) for s = lb
if tau2 < 0; flag_kcase = 1; end % maximize k(t) for s = ub
        
if tau1 < 0 % hard case 2 can only occur towards lb
  if vaaup < zerotol % We have the 'almost' hard case
    % try to take advantage of hard case, case 2 to solve for opt explicity using only pcg
    flag_shiftA = -1;  % indicates a shift by lambdaup*B for Ahandle only
    flag_posshc2 = 1;    % indicates possible hard case, case 2 
    lambda = lambdaup;

       %% save values for the deflation if multiple min eigenvalue
    alphall = []; vall = [];  lambdall = [];       % save alpha, va, lambda
    nAn = nA/sqrt(n); % nA=||A||_1, the value used for the deflation

    while flag_posshc2,  
      % while we still have the possible hard case (case 2)
      if vaaup < zerotol,   % indicates possible hard case, case 2
        alphall = [alphall; nAn];  % save the value of the alpha
        vall = [vall, vaup];   % save the value of va, eigenvector, ensure norm 1
        lambdall = [lambdall, lambda]; % save the value of lambda
        opts.initialvec = randn(n,1);
        opts.size = n;
        opts.maxit = 10000;
        eigtime = tic;
        [lambda,vaup] = eigifp('Ahandle',1,opts); % Compute the smallest eigenvalue after shift and deflation
        fprintf(' eig diff = %g, time = %g\n', norm(Ahandle(vaup) - lambda*vaup), toc(eigtime))
        vaup = vaup/norm(vaup);     % scaling

        vaaup = vaup'*a/na;
        if vaaup < 0 % loop to ensure the inner product is positive
          vaup = -vaup;
          vaaup = -vaaup;
        end 
             
      else
        flag_posshc2 = 0;            
      end 
      
      if lambda >= zerotol
        break
      end
    end

    if flag_posshc2   % possible hard case, case 2 ... check norm of x found using pcg
      disp('**POSSIBLE** hard case, case 2; using pcg to find x');

      [xtemp, dummy] = pcg('Ahandle',a,tol,1000);  % solve the x(\lambda) 
      Bxtemp = Bhandle(xtemp);
      tau1 = xtemp'*Bxtemp - lb; % compute tau - step length to boundary              

      if tau1 <= 0,                  % we do have the hard case
	      disp('HARD CASE, CASE 2; taking step to lower boundary');
        vb = vall(:,1)/sqrt(vall(:,1)'*Bhandle(vall(:,1))); % Since vall(:,1)) is in the nullspace of A - lambdaup*B, 
                                                      % A is pd and lambdaup > 0, it follows that vall(:,1)'*B*vall(:,1) > 0
        b = vb'*Bxtemp;
        c = tau1; % The scaling in vb gives a quadratic equation with a = 1
        if b < 0
          stepsize = -b + sqrt(b^2 - c);
        else
          stepsize = c/(-b - sqrt(b^2 - c));
        end
        xopt = xtemp + stepsize*vb;
        lopt = lambdaup;
        flag_shiftA = 0;
        fopt = xopt'*Ahandle(xopt) - 2*a'*xopt;
              
        flag_case = 1;
        output.y = xopt;
        output.fval = fopt;
        output.lopt = lopt;
        constraint = xopt'*Bhandle(xopt);
        output.feas = constraint - max(min(constraint,ub),lb);
        s = lb;
        return 
      else                   
        flag_posshc2 = 0;
      end
               
    else  % x(\lambda) > lb, not the hard case (2)
      flag_posshc2 = 0;
    end
  end        
end

if tau2 < 0 
  if vaalow < zerotol % We have the 'almost' hard case
    % try to take advantage of hard case, case 2 to solve for opt explicity using only pcg
    flag_shiftA = 1;  % indicates a shift by lambdalow*B for Ahandle only
    flag_posshc2 = 1;    % indicates possible hard case, case 2  
    lambda = lambdalow;

       %% save values for the deflation if multiple min eigenvalue
    alphall = []; vall = [];  lambdall = [];       % save alpha, va, lambda
    nAn = nA/sqrt(n); % nA=||A||_1, the value used for the deflation

    while flag_posshc2,  
      % while we still have the possible hard case (case 2)
      if vaalow < zerotol,   % indicates possible hard case, case 2
        alphall = [alphall; nAn];  % save the value of the alpha
        vall = [vall, valow];   % save the value of va, eigenvector, ensure norm 1
        lambdall = [lambdall, lambda]; % save the value of lambda
        opts.initialvec = randn(n,1);
        opts.size = n;
        opts.maxit = 10000;
        eigtime = tic;
        [lambda,valow] = eigifp('Ahandle',1,opts); % Compute the smallest eigenvalue after shift and deflation
        fprintf(' eig diff = %g, time = %g\n', norm(Ahandle(valow) - lambda*valow), toc(eigtime))
        valow = valow/norm(valow);     % scaling

        vaalow = valow'*a/na;
        if vaalow < 0 % loop to ensure the inner product is positive
          valow = -valow;
          vaalow = -vaalow;
        end 
             
      else
        flag_posshc2 = 0;            
      end 
      
      if lambda >= zerotol
        break
      end
    end

    if flag_posshc2   % possible hard case, case 2 ... check norm of x found using pcg
      disp('**POSSIBLE** hard case, case 2; using pcg to find x');

      [xtemp, dummy] = pcg('Ahandle',a,tol,1000);  % solve the x(\lambda) 
      Bxtemp = Bhandle(xtemp);
      tau2 = ub - xtemp'*Bxtemp; % compute tau - step length to boundary              

      if tau2 <= 0,                  % we do have the hard case
	      disp('HARD CASE, CASE 2; taking step to upper boundary');
        vb = vall(:,1)/sqrt(-(vall(:,1)'*Bhandle(vall(:,1))));  % Since vall(:,1)) is in the nullspace of A - lambdalow*B, 
                                                      % A is pd and lambdalow < 0, it follows that vall(:,1)'*B*vall(:,1) < 0
        b = -vb'*Bxtemp;
        c = tau2; % scaling of vb was done so that a = -1; rearrange the quadratic equation to have leading coefficient 1
        if b < 0
          stepsize = -b + sqrt(b^2 - c);
        else
          stepsize = c/(-b - sqrt(b^2 - c));
        end
        xopt = xtemp + stepsize*vb;
        lopt = lambdalow;
        flag_shiftA = 0;
        fopt = xopt'*Ahandle(xopt) - 2*a'*xopt;
              
        flag_case = 1;
        output.y = xopt;
        output.fval = fopt;
        output.lopt = lopt;
        constraint = xopt'*Bhandle(xopt);
        output.feas = constraint - max(min(constraint,ub),lb);
        s = ub;
        return 
      else                   
        flag_posshc2 = 0;
      end
               
    else  % x(\lambda) < ub, not the hard case (2)
      flag_posshc2 = 0;
    end
  end      
end
    

flag_case = 0;
output = 0;
flag_shiftA = 0; % forget the shift, unless the generalized eigenspaces are orthogonal. 
                 
% maximize k(t) with s = ub or lb, according to flag_kcase
if flag_kcase == 0
  s = lb;
else
  s = ub;
end

end

function [v,v_init,dummy,normCv] = GTRS_ktsubproblem

% Solves
%
% min x'*C*x s.t. x'*D*x = 1
%
% Assumes C > 0.

[dummy,v] = eigifp('Dhandle','Chandle',1,opts);
v_init = v;
% normCv = sqrt(v'*Chandle(v));
% v = v/normCv;

v = v/sqrt(-dummy);

end

function tmid = trig_interpol(gindex,bindex,ptgood,ptbad,high,low,tmid)

if bindex > 0 && gindex > 0
  
  % k'(t) from the last point on either side
  tbslope = ptbad(1,3);
  tgslope = ptgood(1,3); 
  
  % k(t) from the last point on either side
  ktb = ptbad(1,4);
  ktg = ptgood(1,4);
  
  ttempt = (tbslope*ptbad(1,1) + ktg - tgslope*ptgood(1,1) - ktb)/(tbslope - tgslope);
  
  if low <= ttempt && ttempt <= high
    tmid = ttempt;
%     fprintf('triang used\n')
  else
%     fprintf('triang fails: outside\n')
  end
end

end

function [tmid,high,low] = vertical_cut(gindex,bindex,ptgood,ptbad,high,low,tmid)

if gindex > 0 && bindex > 0
  
  ktg = ptgood(1,4);
  ktb = ptbad(1,4);
  
  if ktg > ktb % value on good side > bad side
    tmplow = ptbad(1,1) + (ktg - ktb)/ptbad(1,3);
    
    if low <= tmplow && tmplow <= high
      low = tmplow;
%       fprintf(' vertical cut: interval reduced\n')
      if tmid < low
        tmid = low;
%         fprintf(' vertical cut used\n')
      end
    end
  else
    tmphigh = ptgood(1,1) + (ktb - ktg)/ptgood(1,3);
    
    if low <= tmphigh && tmphigh <= high
      high = tmphigh;
%       fprintf(' vertical cut: interval reduced\n')
      if tmid > high
        tmid = high;
%         fprintf(' vertical cut used\n')
      end
    end
  end
end

end

function tmid = inv_interpol(gindex,bindex,ptgood,ptbad,side,ss1,high,low,tmid)

if side < 0 && bindex > 1 % inverse interpolation
  tint = ptbad(1:2,1); 
  vint = ptbad(1:2,2);
  
  f1 = -1/vint(1) + ss1; 
  f2 = -1/vint(2) + ss1;
    
  slope = (tint(1) - tint(2))/(f1 - f2);
  ttemp = tint(1) - f1*slope; % linear inverse interpolation
  if (low <= ttemp) && (ttemp <= high)
    tmid = ttemp;
%     fprintf(' inv interpol used\n')
  else
%     fprintf(' inv interpol failed\n')
  end
else
  if side > 0 && gindex > 1 % inverse interpolation
    tint = ptgood(1:2,1); 
    vint = ptgood(1:2,2);
    
    f1 = -1/vint(1) + ss1; 
    f2 = -1/vint(2) + ss1;
  
    slope = (tint(1) - tint(2))/(f1 - f2);
    ttemp = tint(1) - f1*slope; % linear inverse interpolation
    if (low <= ttemp) && (ttemp <= high)
      tmid = ttemp;
%       fprintf(' inv interpol used\n')
    else
%       fprintf(' inv interpol failed\n')
    end
  end
end

end

end