classdef liu_tests < meq_test
  % Tests for LIUQE
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    tok = 'ana';
    t = 0;
    verbose = 0;
    tol = 1e-5;
    L,LX;
  end

  properties(ClassSetupParameter)
    config = {'limited','diverted','secondaryX'};
    iters = struct('off_0',0,'on_48',48);
    itert={0 1};
  end
  
  properties(TestParameter)
    debug={0,1,2,3};
  end
  
  methods(TestClassSetup)
    function setup_LLX(testCase,config,iters,itert)
      switch config
        case 'limited'
          shot = 1;
        case 'diverted'
          shot = 2;
        case 'secondaryX'
          shot = 3;
        otherwise
          error('undefined config %s\n',config)
      end
      tok = testCase.tok; t=testCase.t;
      
      %FBT
      [L,LX] = liu(tok,shot,t,'bfct',@bfabmex,'bfp',[1,2],'iters',iters,...
        'itert',itert,'iterq',itert);
      
      testCase.LX = LX;
      testCase.L = L;
    end
  end
  
  methods(Test,TestTags={'Unit'})
    function test_ag_constraint_exact(testCase)
      % Tests for constraining ag values.
      
      L = testCase.L; %#ok<*PROP,*PROPLC>
      LX = testCase.LX;
      nt = numel(LX.t);
      
      %% Initial LIUQE run with bfp = [1,2]
      LY0 = liut(L,LX,'debug',testCase.verbose);

      %% New LIUQE run with bfp = [1,3] and constraining last basis function coefficient to 0
      P = L.P;
      G = L.G;
      P.bfp = [1,3];
      P.wag = [0;0;0;Inf]; % This is enough to enable the exact constraint
      L = liuc(P,G);
      LX.ag = 0*ones(L.ng,nt);

      LY1 = liut(L,LX,'debug',testCase.verbose);
      
      % Verifications
      mask = isinf(P.wag);
      tol = testCase.tol*max(abs(LY0.ag));
      testCase.verifyEqual(LY1.ag(mask),LX.ag(mask),'Constrained ag value does not match its input value');
      testCase.verifyEqual(LY1.ag(1:3) ,LY0.ag(1:3),'AbsTol',tol,'Results of two LIUQE runs do not match');
     
    end
    
    function test_ag_constraint_approx(testCase)
      % Tests for constraining ag values.
      
      L = testCase.L;
      LX = testCase.LX;
      nt = numel(LX.t);

      %% Runs with large increasing weights for last basis function coefficient
      P = L.P;
      G = L.G;
      P.bfp = [1,3];
      P.wag = [0;0;0;1e-2]; % This is enough to dominate the cost function
      L = liuc(P,G);
      LX.ag = 0*ones(L.ng,nt);

      LY1 = liut(L,LX,'debug',testCase.verbose);
      
      P.wag = P.wag*100;
      L = liuc(P,G);

      LY2 = liut(L,LX,'debug',testCase.verbose);
      
      % Verifications
      mask = logical(P.wag);
      testCase.verifyLessThan(abs(LY2.ag(mask) - LX.ag(mask)),abs(LY1.ag(mask) - LX.ag(mask)),...
        'Constrained ag value does not decrease with large increasing weights');
     
    end
    
    function debug_test(testCase,debug)
      % check that all debug levels work
      L = testCase.L; LX = testCase.LX;
      if L.P.shot == 1
        L.P.debug = debug;
        L.P.itera = 2; % shorten
        LY = liut(L,LX);
        testCase.verifyTrue(true);
      end
    end
  end
end
