classdef rzp_fast_test < meq_test
  %% Test class for rzpfastA
  % Test the working flow of rzpfastA with different input parameters,
  % its mexification and compare the eigenvalues calculated with this fast
  % version versus those calculated with the slow one.
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0; % 1 to make verbose output, 2 to make plots

    selu    = 'e';
    selx    = 'X';    
    insrc   = 'fbt';
    iterq   = 20;
    agcon   = {'Ip','bp','li'};
    rtol    = 2.e-2; %tolerance on the elements relative difference 
    n_eigen = 5;     %number of eigenvalues for the fast-slow comparison
    
    tok,shot,t,nu,cde,Lliu,LXliu;
  end
  
  properties (ClassSetupParameter)
    tok_            = {'ana','tcv'};  %tokamak type
    tnum            = {1,10};         %number of time points
    shot_time_index = {1,2,3};        %index for the shot number (different lists are used in anamak and tcv case)
  end
  
  properties (MethodSetupParameter)
    cde_type  = struct('nu_max', ''   , 'no_cde', ''  , 'cde_OhmTorRigid', 'OhmTor_rigid'); 
    vessel_nu = struct('nu_max', 'max', 'no_cde', 30  , 'cde_OhmTorRigid',   30          );
  end
  
  methods (Static)
    function A_fast = fast_compacted_input(L,LX)
      A_fast = rzpfastA(LX.Ia,LX.Iu,LX.Iy,LX.Ip,LX.rIp,LX.bp,...
        L.Mey,L.Mee,L.Re,LX.Rp,L.RBrye,L.RBzye,L.drx,L.dzx,L.ny,L.np,L.ne,L.icde);
    end
  end
  
  methods (TestClassSetup)
    function setup_parameters(testCase,tok_,tnum,shot_time_index)
      testCase.tok  = tok_;
      if strcmp(tok_,'ana')  
        shot_list     = [1,2,4]; 
        time_list     = {linspace(0.005,0.015,tnum),linspace(0.005,0.015,tnum),linspace(0.005,0.015,tnum)};
        testCase.shot = shot_list(shot_time_index);
        testCase.t    = time_list{shot_time_index};
      else 
        shot_list     = [61400,69393,63783]; 
        time_list     = {linspace(0.99,1,tnum),linspace(0.99,1,tnum),linspace(0.99,1,tnum)};
        testCase.shot = shot_list(shot_time_index);
        testCase.t    = time_list{shot_time_index};
      end

      [testCase.Lliu,testCase.LXliu] = fgs(testCase.tok,testCase.shot,testCase.t,'selu',testCase.selu,'selx',testCase.selx,'insrc',testCase.insrc);
    end
  end
  
  methods (TestMethodSetup, ParameterCombination='sequential')
    function setup_nu(testCase,vessel_nu,cde_type)
      testCase.cde = cde_type;
      if strcmp(vessel_nu,'max')
        if strcmp(testCase.tok,'ana')
          testCase.nu = 200;
        else
          testCase.nu = 256;
        end
      else
        testCase.nu = vessel_nu;
      end
    end
  end
  
  methods (Test, TestTags = {'rzp'}, ParameterCombination='sequential')   
    function code_test_parameter(testCase)    
      %% Test that rzpfastA is operating
      case_string = sprintf('fast A calculation with %s_%u_%.5g',testCase.tok,testCase.shot,testCase.t(1));
      P_args = {'insrc',testCase.insrc,'selu',testCase.selu,'selx',testCase.selx,...
        'agcon',testCase.agcon,'debug',testCase.verbosity,'cde',testCase.cde,'nu',testCase.nu};
      P_error = sprintf('%s: %s, %s: %s, %s: %d',P_args{1:2},P_args{11:14});  
      
      L = rzp(testCase.tok,testCase.shot,[],P_args{:});
      % Get inputs data for simulation
      LX = meqxconvert(testCase.Lliu,testCase.LXliu,L);
      LX = rzpx(testCase.tok,testCase.t,L,LX);
      % compute linearization and initial condition
      L = rzpl(L,meqxk(LX,1));
      
      try
        for ii = 1:numel(testCase.t)
         testCase.fast_compacted_input(L,meqxk(LX,ii));   
        end
      catch ME
        fprintf('\n %s error caused in %s in rzpfastA by parameters %s \n',...
          ME.message,case_string,P_error);
          rethrow(ME);
      end 
    end
    
    function mexify_test(testCase)    
      %% Test that rzpfastAmexify is operating    
      if numel(testCase.t) == 1 && testCase.nu == 30 && strcmp(testCase.cde,'')
        case_string = sprintf('fast A calculation with %s_%u_%.5g',testCase.tok,testCase.shot,testCase.t(1));
        P_args = {'insrc',testCase.insrc,'selu',testCase.selu,'selx',testCase.selx,...
          'agcon',testCase.agcon,'debug',testCase.verbosity,'cde',testCase.cde,'nu',testCase.nu};
        P_error = sprintf('%s: %s, %s: %s, %s: %d',P_args{1:2},P_args{11:14});  

        [L] = rzp(testCase.tok,testCase.shot,testCase.t,P_args{:});
        store_dir = sprintf('%s/meq/codegen',tempname);
        out_dir   = store_dir;
        out_file  = fullfile(out_dir,['rzpfastA_mex.',mexext]);
        try
          if exist(out_dir,'dir')
            rmdir(out_dir,'s');
          end
          rzpfastAmexify(L,store_dir,out_dir);  
        catch ME
          fprintf('\n %s error caused in %s in rzpfastAmexify by parameters %s \n',...
            ME.message,case_string,P_error);
          rethrow(ME);
        end
        rmdir(out_dir,'s');
      end
    end
    
    function eigen_test(testCase)
      %% Test fast vs slow eigenvalue calculation
      if testCase.nu == 30 && numel(testCase.t) == 1
        case_string = sprintf('fast A calculation with %s_%u_%.5g',testCase.tok,testCase.shot,testCase.t(1));
        P_args = {'iterq',testCase.iterq,'insrc',testCase.insrc,'selu',testCase.selu,'selx',testCase.selx,...
          'agcon',testCase.agcon,'debug',testCase.verbosity,'cde',testCase.cde,'nu',testCase.nu};
        P_error = sprintf('%s: %s, %s: %s, %s: %d',P_args{1:2},P_args{11:14}); 
      
        L = rzp(testCase.tok,testCase.shot,[],P_args{:});
        % Get inputs data for simulation
        LX = meqxconvert(testCase.Lliu,testCase.LXliu,L);
        LX = rzpx(testCase.tok,testCase.t,L,LX);
        % compute linearization and initial condition
        L = rzpl(L,meqxk(LX,1));

        eig_fast = esort(eig(testCase.fast_compacted_input(L,LX)));
        
        sys = fgess(L); 
        eig_slow = esort(eig(sys.A));
        
        if eig_slow(1) > 0 
          testCase.verifyEqual( eig_fast(1:testCase.n_eigen), eig_slow(1:testCase.n_eigen),...
            'RelTol', testCase.rtol,sprintf('Eigen fast vs slow tolerance exceeded in %s by paramters %s',case_string,P_error));
        end
      end
    end
  end
end
