
% for pdA and indB

clc
clear

density = 1e-2;
repeats = 10;
% 0 for easy case, 1 for hard case 1, 2 for hard case 2, 3 for interior solution

% options for eigifp
opts.useprecon = 'NO';
opts.disp = 0;
opts.maxit = 10000;
fprintf(' Generating\n')

fid1 = fopen('indBpdAe.txt','w');
fid3 = fopen('indBpdAh2.txt','w');

for n = 1000*[10 15 20]
  
  for cases = [0 2]
    
    iter1 = 0; iter3 = 0;
    time1 = 0; time3 = 0;
    fval1 = 0; fval3 = 0;
    feas1 = 0; feas3 = 0;
    
    repeat = 1;
    while repeat <= repeats
      
      tst = tic;
      A = sprandsym(n,density,0.1,2); % very low density
      B = sprandsym(n,density);
      
      opts.initialvec = sum(A)';
      mu = eigifp(-B,A,1,opts);
      lambdaup = -1/mu;
      mu = eigifp(B,A,1,opts);
      lambdalow = 1/mu;
      fprintf(' Done generation, time %g\n', toc(tst))
      
      A1 = A;
      B1 = B;
      fprintf(' Done generation, time %g\n', toc(tst))
      
      if cases == 0  
        
        lambda = (lambdalow + lambdaup)/2; % avoid hard case. 
        x = randn(n,1)/10;   
        a = (A1 - lambda*B1)*x;
        s = x'*B1*x;
        % scale B to get t0 and lambda0
        B = B1/(-abs(s));
        A = A1;
        lambda0 = lambda*(-abs(s));
        t0 = 2*a'*x - x'*A*x;
        s = randn(1)*sqrt(n);
        lb = s;
        ub = s + 1;
      
      elseif cases == 2
        
        tmp = rand(1);
        if tmp > 0.5 % lower boundary hard case 2
          lambda = lambdalow;
          x = randn(n,1)/10;
          a = (A1 - lambda*B1)*x;
          s = x'*B1*x;
          % scale B to get t0 and lambda0
          B = B1/(-abs(s));
          A = A1;
          lambda0 = lambda*(-abs(s));
          t0 = 2*a'*x - x'*A*x;
          s = (rand(1) + 1)*sqrt(n);
          lb = s;
          ub = s + 1;
        else
          lambda = lambdaup; % upper boundary hard case 2
          x = randn(n,1)/10;
          a = (A1 - lambda*B1)*x;
          s = x'*B1*x;
          % scale B to get t0 and lambda0
          B = B1/(-abs(s));
          A = A1;
          lambda0 = lambda*(-abs(s));
          t0 = 2*a'*x - x'*A*x;
          s = -(rand(1) + 1)*sqrt(n);
          lb = s;
          ub = s + 1;
        end
        
      end
      
      fprintf(' RW \n')
      [y,iter,fval,feas,dval,k_slope,time,lambda,side] = gtrsalgo_ul_indB(A,a,B,lb,ub,1e-13,t0,lambda0);
      feas = y'*B*y - side;
    
      if side == side
        iter1 = iter1 + iter;
        time1 = time1 + time;
        fval1 = fval1 + fval;
        feas1 = feas1 + abs(feas);
        
        fprintf(' Newton \n')
        [u,lambda,fval,feas,iter,time] = Newton_gtrs_indB(A,a,B,side,1e-12,10);
        iter3 = iter3 + iter;
        time3 = time3 + time;
        fval3 = fval3 + fval;
        feas3 = feas3 + abs(feas);
        repeat = repeat + 1;
      else % cancel the instance and regenerate
        fprintf('\n\n REGENERATE !!!\n\n')
      end
      
    end
  
    if cases == 0; 
      fprintf(fid1,'%5d & %3.0f/%6.2f/%10.9e/%2.1e &', n, iter1/repeats, time1/repeats, fval1/repeats, feas1/repeats);
      fprintf(fid1,'%3.0f/%6.2f/%10.9e/%2.1e\n', iter3/repeats, time3/repeats, fval3/repeats, feas3/repeats);
    elseif cases == 2;
      fprintf(fid3,'%5d & %3.0f/%6.2f/%10.9e/%2.1e &', n, iter1/repeats, time1/repeats, fval1/repeats, feas1/repeats);
      fprintf(fid3,'%3.0f/%6.2f/%10.9e/%2.1e\n', iter3/repeats, time3/repeats, fval3/repeats, feas3/repeats);
    end
    
  end

end

fclose(fid1);
fclose(fid3);