classdef (SharedTestFixtures={mexm_fixture}) ...
 asxy_test < meq_test
  % (See also analytical_map_test.m) 
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    data;
    flux_abs_tol = 1e-2;
    hess_abs_tol = 1e-1;
    grid_rel_tol = 1;
    verbosity = 0;
  end
  
  properties (ClassSetupParameter)
    type = {'SND',...
            'DND',...
            'Lim',...
            'Lim-X',...           % limited with a non-limiting X-point
            'Boundary-X',...      % X-point close to the computational boundary
            'Double-Snowflake-Minus',...
            'Double-Snowflake-Plus',...
            'Doublet',...
            'Droplet',...
            'Doublet-div-nomantle',...
            'Doublet-div',...
            'Triplet',...
            'Triplet-madness',...
           };
    gridsize = {[31,32],  [61,64]};
  end
  
  methods(TestClassSetup)
    function asxyTestSetup(testCase,type,gridsize)
      S = testCase.generate_flux_map(type,1);
      
      L = fgs('ana',1,[],'nr',gridsize(1)-1,'nz',gridsize(2),...
              'ri',0.55,'ro',1.45,'zl',-0.9,'zu',0.9,...
              'cappav',2,'icsint',true); 
      
      Fx = S.Fxh(L.rrx,L.zzx);
      
      testCase.data.S  = S;
      testCase.data.Fx = Fx;
      testCase.data.L  = L;
    end
  end
  
  properties(TestParameter)
    alternate_orientation = {'ud','ud+lr','lr','T','T+ud','T+ud+lr','T+lr'};
    method = {'asxymex','asxymexm'};
    refining = struct('off',false,'on',true);
  end
  
  methods(Test, TestTags = {'Unit'})
    function test_asxy(testCase,method,refining)
      % Compare results with fminsearch for analytical flux maps.
      S  = testCase.data.S;
      Fx = testCase.data.Fx;
      L  = testCase.data.L;
      
      %% Tests for asxymex - Localise extrema and saddle points of a 2D-map
      % [XA,YA,FA,DX2FA,DY2FA,DXYFA,XS,YS,FS,DX2FS,DY2FS,DXYFS] = ASXYMEX(F,X,Y,DX,DY,IDX,IDY,L)
      
      % Localise extrema and saddle points of the Fx map
      [zA,rA,FA,dz2FA,dr2FA,drzFA,zX,rX,FX,dz2FX,dr2FX,drzFX] = feval(method,Fx,L.G.zx,L.G.rx,2,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx);

      if refining
        [zA,rA,FA,dz2FA,dr2FA,drzFA,zX,rX,FX,dz2FX,dr2FX,drzFX] = asxycs(Fx,zA,rA,zX,rX,L);
      end
      
      % Verify it found the correct number of axes and X-points
      
      testCase.assertEqual(numel(zA),S.nA,'Numbers of axes do not match');
      testCase.assertEqual(numel(zX),S.nX,'Numbers of X-points do not match');
      
      % Compare with fminsearch results
      [zA_,rA_,FA_,dz2FA_,dr2FA_,drzFA_] = deal(NaN(S.nA,1));
      for iA = 1:S.nA
        [x] = fminsearch(@(x) S.drFxh(x(1),x(2)).^2 + S.dzFxh(x(1),x(2)).^2,[rA(iA),zA(iA)]);
        zA_(iA) = x(2);
        rA_(iA) = x(1);
        FA_(iA) = S.Fxh(x(1),x(2));
        dz2FA_(iA) = S.dz2Fxh(x(1),x(2));
        dr2FA_(iA) = S.dr2Fxh(x(1),x(2));
        drzFA_(iA) = S.drzFxh(x(1),x(2));
      end
      [zX_,rX_,FX_,dz2FX_,dr2FX_,drzFX_] = deal(NaN(S.nX,1));
      for iX = 1:S.nX
        [x] = fminsearch(@(x) S.drFxh(x(1),x(2)).^2 + S.dzFxh(x(1),x(2)).^2,[rX(iX),zX(iX)]);
        zX_(iX) = x(2);
        rX_(iX) = x(1);
        FX_(iX) = S.Fxh(x(1),x(2));
        dz2FX_(iX) = S.dz2Fxh(x(1),x(2));
        dr2FX_(iX) = S.dr2Fxh(x(1),x(2));
        drzFX_(iX) = S.drzFxh(x(1),x(2));
      end

      if refining
        tolf = 2e-1;
      else
        tolf = 1e+0;
      end
      
      % Axis checks
      testCase.verifyEqual(   rA,   rA_, 'AbsTol',L.drx*testCase.grid_rel_tol*tolf,   'rA result failure');
      testCase.verifyEqual(   zA,   zA_, 'AbsTol',L.dzx*testCase.grid_rel_tol*tolf,   'zA result failure');
      testCase.verifyEqual(   FA,   FA_, 'AbsTol',      testCase.flux_abs_tol*tolf,   'FA result failure');
      testCase.verifyEqual(dr2FA,dr2FA_, 'AbsTol',      testCase.hess_abs_tol*tolf,'dr2FA result failure');
      testCase.verifyEqual(dz2FA,dz2FA_, 'AbsTol',      testCase.hess_abs_tol*tolf,'dz2FA result failure');
      testCase.verifyEqual(drzFA,drzFA_, 'AbsTol',      testCase.hess_abs_tol*tolf,'drzFA result failure');
      
      % X-point checks
      testCase.verifyEqual(   rX,   rX_, 'AbsTol',L.drx*testCase.grid_rel_tol*tolf,   'rX result failure');
      testCase.verifyEqual(   zX,   zX_, 'AbsTol',L.dzx*testCase.grid_rel_tol*tolf,   'zX result failure');
      testCase.verifyEqual(   FX,   FX_, 'AbsTol',      testCase.flux_abs_tol*tolf,   'FX result failure');
      testCase.verifyEqual(dr2FX,dr2FX_, 'AbsTol',      testCase.hess_abs_tol*tolf,'dr2FX result failure');
      testCase.verifyEqual(dz2FX,dz2FX_, 'AbsTol',      testCase.hess_abs_tol*tolf,'dz2FX result failure');
      testCase.verifyEqual(drzFX,drzFX_, 'AbsTol',      testCase.hess_abs_tol*tolf,'drzFX result failure');
      
      if testCase.verbosity
        figure;
        contour(L.G.rx,L.G.zx,Fx,50,'-r')
        hold on, plot(rX_,zX_,'xk',rX,zX,'sg'), hold off, axis equal
        txt = sprintf('%dx%d',numel(L.G.rx),numel(L.G.zx));
        title(['X point - ',txt])
        drawnow
        shg
      end
      
    end
    
    function test_asxy_orientation(testCase,alternate_orientation,method)
      % Check that the result is unchanged if the orientation is changed
      Fx = testCase.data.Fx;
      L  = testCase.data.L;
      
      %% Tests for asxymex - Localise extrema and saddle points of a 2D-map
      % [XA,YA,FA,DX2FA,DY2FA,DXYFA,XS,YS,FS,DX2FS,DY2FS,DXYFS] = ASXYMEX(F,X,Y,DX,DY,IDX,IDY,L)
      
      % Localise extrema and saddle points of the Fx map
      [zA_,rA_,~,~,~,~,zX_,rX_,~,~,~,~] = feval(method,Fx,L.G.zx,L.G.rx,2,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx);
      
      switch alternate_orientation
        case 'ud'
          [zA,rA,~,~,~,~,zX,rX,~,~,~,~] = feval(method,flipud(Fx),L.G.zx(end:-1:1),L.G.rx,2,-L.dzx,L.drx,-L.idzx,L.idrx,flipud(L.Oasx));
        case 'ud+lr'
          [zA,rA,~,~,~,~,zX,rX,~,~,~,~] = feval(method,rot90(Fx,2),L.G.zx(end:-1:1),L.G.rx(end:-1:1),2,-L.dzx,-L.drx,-L.idzx,-L.idrx,rot90(L.Oasx,2));
        case 'lr'          
          [zA,rA,~,~,~,~,zX,rX,~,~,~,~] = feval(method,fliplr(Fx),L.G.zx,L.G.rx(end:-1:1),2,L.dzx,-L.drx,L.idzx,-L.idrx,fliplr(L.Oasx));
        case 'T'
          [rA,zA,~,~,~,~,rX,zX,~,~,~,~] = feval(method,Fx.',L.G.rx,L.G.zx,2,L.drx,L.dzx,L.idrx,L.idzx,L.Oasx.');
        case 'T+ud'
          [rA,zA,~,~,~,~,rX,zX,~,~,~,~] = feval(method,flipud(Fx.'),L.G.rx(end:-1:1),L.G.zx,2,-L.drx,L.dzx,-L.idrx,L.idzx,flipud(L.Oasx.'));
        case 'T+ud+lr'
          [rA,zA,~,~,~,~,rX,zX,~,~,~,~] = feval(method,rot90(Fx.',2),L.G.rx(end:-1:1),L.G.zx(end:-1:1),2,-L.drx,-L.dzx,-L.idrx,-L.idzx,rot90(L.Oasx.',2));
        case 'T+lr'
          [rA,zA,~,~,~,~,rX,zX,~,~,~,~] = feval(method,fliplr(Fx.'),L.G.rx,L.G.zx(end:-1:1),2,L.drx,-L.dzx,L.idrx,-L.idzx,fliplr(L.Oasx.'));
      end
      testCase.verifyEqual(numel(zA),numel(zA_),sprintf('Numbers of axes do not match for orientation %s',alternate_orientation));
      testCase.verifyEqual(numel(zX),numel(zX_),sprintf('Numbers of X-points do not match for orientation %s',alternate_orientation));
      % Order might have changed
      iA = testCase.closest_pairs(rA_,zA_,rA,zA);
      rA = rA(iA);zA = zA(iA);
      iX = testCase.closest_pairs(rX_,zX_,rX,zX);
      rX = rX(iX);zX = zX(iX);
      %
      testCase.verifyEqual(rA,rA_, 'AbsTol',L.drx*testCase.grid_rel_tol,sprintf('rA result failure for orientation %s',alternate_orientation));
      testCase.verifyEqual(zA,zA_, 'AbsTol',L.drx*testCase.grid_rel_tol,sprintf('zA result failure for orientation %s',alternate_orientation));
      testCase.verifyEqual(rX,rX_, 'AbsTol',L.drx*testCase.grid_rel_tol,sprintf('rX result failure for orientation %s',alternate_orientation));
      testCase.verifyEqual(zX,zX_, 'AbsTol',L.drx*testCase.grid_rel_tol,sprintf('zX result failure for orientation %s',alternate_orientation));
      
    end
    
    function test_asxy_implementation(testCase)
      % Compare C and MATLAB implementation
      
      Fx = testCase.data.Fx;
      L  = testCase.data.L;
      
      %% Tests for asxymex - Localise extrema and saddle points of a 2D-map
      % [XA,YA,FA,DX2FA,DY2FA,DXYFA,XS,YS,FS,DX2FS,DY2FS,DXYFS] = ASXYMEX(F,X,Y,DX,DY,IDX,IDY,L)
      
      % Localise extrema and saddle points of the Fx map - using mex
      [zA_,rA_,FA_,dz2FA_,dr2FA_,drzFA_,zX_,rX_,FX_,dz2FX_,dr2FX_,drzFX_] = asxymex(Fx,L.G.zx,L.G.rx,2,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx);
      
      % Localise extrema and saddle points of the Fx map - using mexm
      [zA ,rA ,FA ,dz2FA ,dr2FA ,drzFA ,zX ,rX ,FX ,dz2FX ,dr2FX ,drzFX ] = asxymexm(Fx,L.G.zx,L.G.rx,2,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx);
      
      testCase.verifyEqual(numel(zA),numel(zA_),'Numbers of axes do not match');
      testCase.verifyEqual(numel(zX),numel(zX_),'Numbers of X-points do not match');
      
      % Axis checks
      testCase.verifyEqual(   rA,   rA_, 'AbsTol',L.drx*testCase.grid_rel_tol,   'rA result failure');
      testCase.verifyEqual(   zA,   zA_, 'AbsTol',L.dzx*testCase.grid_rel_tol,   'zA result failure');
      testCase.verifyEqual(   FA,   FA_, 'AbsTol',      testCase.flux_abs_tol,   'FA result failure');
      testCase.verifyEqual(dr2FA,dr2FA_, 'AbsTol',      testCase.hess_abs_tol,'dr2FA result failure');
      testCase.verifyEqual(dz2FA,dz2FA_, 'AbsTol',      testCase.hess_abs_tol,'dz2FA result failure');
      testCase.verifyEqual(drzFA,drzFA_, 'AbsTol',      testCase.hess_abs_tol,'drzFA result failure');
      
      % X-point checks
      testCase.verifyEqual(   rX,   rX_, 'AbsTol',L.drx*testCase.grid_rel_tol,   'rX result failure');
      testCase.verifyEqual(   zX,   zX_, 'AbsTol',L.dzx*testCase.grid_rel_tol,   'zX result failure');
      testCase.verifyEqual(   FX,   FX_, 'AbsTol',      testCase.flux_abs_tol,   'FX result failure');
      testCase.verifyEqual(dr2FX,dr2FX_, 'AbsTol',      testCase.hess_abs_tol,'dr2FX result failure');
      testCase.verifyEqual(dz2FX,dz2FX_, 'AbsTol',      testCase.hess_abs_tol,'dz2FX result failure');
      testCase.verifyEqual(drzFX,drzFX_, 'AbsTol',      testCase.hess_abs_tol,'drzFX result failure');
      
    end
  end
  
  
  methods (Static)
    function i2 = closest_pairs(r1,z1,r2,z2)
      i2 = 1:numel(r2);
      for i1 = 1:numel(r1)-1
        [~,ii] = min((r1(i1)-r2(i2(i1:end))).^2+(z1(i1)-z2(i2(i1:end))).^2);
        tmp = i2(i1);
        i2(i1) = i2(i1+ii-1);
        i2(i1+ii-1) = tmp;
      end
    end
  end
  
end
