classdef (SharedTestFixtures={mexm_fixture}) ipm1_test < meq_test
  % tests for IPM1
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    x,H,f,a,UH,xi,zi,Nx
    itertol = 1e-3;
    errtol = 1e-2;
    maxiter = 150;
    verbosity = 0;
  end
  
  properties(TestParameter)
    method = {'ipm1mex','ipm1mexm'};
  end
    
  methods(TestClassSetup)
    function ipm1TestSetup(testCase)
            
      testCase.assumeTrue(~isempty(which('quadprog')), 'quadprog missing, need optimization toolbox for fbtt')
      if ~isempty(which('optimoptions')) % optimoptions not available in Octave
        opt = optimoptions('quadprog','algorithm','interior-point-convex','display','off');
      else
        opt = [];
      end
      
      testCase.Nx = 2;
      xy0 = rand(1,testCase.Nx);
      t = 360*rand;
      ax1 = 1.5 + rand;
      ax2 = 1.5 + rand;
      testCase.a = sign(rand-0.5);
      if not(testCase.a), testCase.a = 1; end
      
      R = [cosd(t) sind(t); -sind(t) cosd(t)];
      HH = R' * diag([1/ax1^2   1/ax2^2]) * R;
      testCase.H = 2 * round((HH+HH')/2, 3);
      testCase.f = round((-[2*xy0(1)/ax1^2    2*xy0(2)/ax2^2] * R)', 3);
      
      testCase.UH = testCase.H(triu(true(size(testCase.H)))); % upper triangular part of H
      testCase.xi = [1;1]*testCase.a;
      testCase.zi = testCase.xi;
      testCase.x  = quadprog(testCase.H,testCase.f,-testCase.a*eye(testCase.Nx),zeros(1,testCase.Nx),[],[],[],[],[],opt);      
    end
  end
  
  methods(Test, TestTags = {'Unit'})
    function ipm1Test(testCase,method)

      [x2,~,ST] = feval(method,testCase.UH,testCase.f,testCase.a,testCase.xi,testCase.zi,testCase.itertol,10);

      testCase.verifyTrue(logical(ST), 'ip1mex procedure returned an error')
      testCase.verifyLessThan(norm(testCase.x-x2),testCase.errtol,'error out of the prescribed bounds, test result failure')
      
      if testCase.verbosity
        clf
        Np = 50;
        xg = testCase.a*linspace(0,5,Np)';
        [xg1,xg2] = meshgrid(xg,xg);
        F = zeros(Np);
        for r = 1:length(xg)
          for c = 1:length(xg)
            xx = [xg1(r,c);xg2(r,c)];
            F(r,c) = xx'*testCase.H*xx/2 + testCase.f'*xx;
          end
        end
        figure
        contour(xg1,xg2,F), grid on, hold on
        plot([0 0], testCase.a*[-0.5 5], '-k', testCase.a*[-0.5 5], [0 0], '-k')
        plot(testCase.x(1),testCase.x(2),'xr',x2(1),x2(2),'^b')
      end
    end
  end 
end
