classdef (SharedTestFixtures={meq_fixture}) ...
    fgslin_full_test < matlab.unittest.TestCase
  % Tests of fgslin_full
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties (TestParameter)
    shot = {61400}
    t = {1};
  end
  
  properties
    verbosity = 1;
    tol = 1e-2;
  end
  
  methods (Test, TestTags = {'fgs'})
    function test_fgs_linearization(testCase,shot,t)
      % checks linearization using fgslin_full (finite differences in all
      % directions and full recomputation of converged FGS equilibrium

      [L,LX,LY] = fgs('tcv',shot,t,'tolF',1e-9,...
        'mkryl',50,'usepreconditioner',0,...
        'selu','e','nu',30,...
        'debug',testCase.verbosity);
      [dLY_dIa,dLY_dIu,dLY_dagcon] = fgslin_full(L,LY);
      LX0 = LX;
      %% check w.r.t. Ia perturbation
      
      LX = LX0;
      dIa = 1e-4*randn(size(LX.Ia)).*LX.Ia;
      dIu = 1e-4*randn(size(LX.Iu));
      LX.Ia = LX.Ia + dIa;
      LX.Iu = LX.Iu + dIu;
      LY1 = fgst(L,LX);
      
      checkfields = {'rA','zA','Fx','FA','FB','PQ','TQ','rB','zB','bp','Ip','qA'};
      for ifield=checkfields
        myfield = ifield{:};
        if isscalar(LY1.(myfield)) || iscolumn(LY1.(myfield))
          linval = [dLY_dIa.(myfield),dLY_dIu.(myfield)] * [dIa;dIu];
        elseif ismatrix(LY1.(myfield))
          linval = sum(cat(3,dLY_dIa.(myfield),dLY_dIu.(myfield)).*reshape([dIa;dIu],1,1,[]),3);
        elseif ndims(LY1.(myfield)) == 3
          linval = sum(cat(4,dLY_dIa.(myfield),dLY_dIu.(myfield)).*reshape([dIa;dIu],1,1,1,[]),4);
        else
          error('invalid field size');
        end
        numval = LY1.(myfield)-LY.(myfield);
        testCase.verifyEqual(numval,linval,'AbsTol',abs(LY.(myfield))*testCase.tol,...
          sprintf('residual too high for %s',myfield));
      end
      %%
      
      LX=LX0; % reset
      L.P.debug=0; % silence these repetitive calls
      % check w.r.t. separate agcon perturbation
      for icon = 1:numel(L.P.agcon)
        mycon = L.P.agcon{icon};
        % special fix for Ip, bp fields (perturb IpD, bpD)
        if any(strcmp(mycon,{'Ip','bp'})); mycon=[mycon,'D']; end %#ok<AGROW> % 
        dcon = 1e-6*LX.(mycon);
        LX.(mycon) = LX.(mycon)+dcon;
        LY1 = fgst(L,LX);
        
        % check that linearization w.r.t. this constraint value gives 1
        testCase.verifyEqual(dLY_dagcon.(mycon)(icon),1,'AbsTol',1e-6,sprintf('linearization around %s constraint is not 1',mycon));
        
        for ifield=checkfields
          myfield = ifield{:};
          if isscalar(LY1.(myfield)) || iscolumn(LY1.(myfield))
            linval = dLY_dagcon.(myfield)(:,icon)*dcon;
          elseif ismatrix(LY.(myfield))
            linval = dLY_dagcon.(myfield)(:,:,icon)*dcon;
          else
            error('invalid field size');
          end
          
          testCase.verifyEqual(LY1.(myfield)-LY.(myfield),linval,...
            'AbsTol',abs(LY.(myfield))*testCase.tol,sprintf('Tolerance exceeded for d%s/d%s',myfield,mycon));
        end
      end
    end
  end
end
