classdef (SharedTestFixtures={mexm_fixture}) iata_test < meq_test
  % Test iatamex
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    tol = 1e-6; % numerical tolerance (for single precision!)
    seed = 0; % randomizer seed
  end
  
  properties (TestParameter)
    n = {'1','2','3','4','5','6','7','8'}; % matrix size for iata
    type = {'single','double'}; % datatype for iata
    method = {'iatamex','iatamexm'};
  end
  
  methods (TestClassSetup)
    function setRandomSeed(testCase)
      rng(testCase.seed); % to always get same 'random' numbers;
    end
  end
  
  methods(Test, TestTags = {'Unit'})
    function test_iata(testCase,n,type,method)
      %% Tests for iata(A1,A2) - very fast matrix inversions
      % should return iA = inv(A1'*A1+A2'*A2)
      n = str2double(n);
      
      % case with A1, A2
      A1 = cast(randn(76,n),type); % size [nr x nj]
      A2 = cast(randn(2,n),type);  % size [ni x nj]
      A=A1'*A1+A2'*A2;
      
      iAmex = feval(method,A1,A2);
      testCase.verifyLessThan(100*eps(type), rcond(A), 'Matrix is badly conditioned, please update tests');
      testCase.verifyLessThan(norm(A*iAmex - eye(n)),testCase.tol, 'iatamex(A1,A2) result failure');
      
      %% case without A2
      A=A1'*A1;
      
      iAmex = feval(method,A1,cast(zeros(0,n),type));
      testCase.assumeLessThan(100*eps(type), rcond(A), 'Matrix is badly conditioned, please update tests');
      testCase.verifyLessThan(norm(A*iAmex - eye(n)),testCase.tol, 'iatamex(A) result failure');

      %%
    end
  end
end
