classdef nfsftUnitTests
  %UNITTEST Summary of this class goes here
  %   Detailed explanation goes here
  
  properties
    perform_exhaustive_tests_flag = @matlab_exhaustive_unit_tests_flag@;

    testcases_online_trafo = { ...
      nfsftTestcaseDelegateOnline(32, 10, 'trafo') ...
      nfsftTestcaseDelegateOnline(32, 100, 'trafo') ...
      nfsftTestcaseDelegateOnline(32, 1000, 'trafo') ...
      nfsftTestcaseDelegateOnline(64, 10, 'trafo') ...
      nfsftTestcaseDelegateOnline(64, 100, 'trafo') ...
      nfsftTestcaseDelegateOnline(128, 10, 'trafo') ...
      nfsftTestcaseDelegateOnline(128, 100, 'trafo') ...
      nfsftTestcaseDelegateOnline(256, 10, 'trafo') ...
      nfsftTestcaseDelegateOnline(256, 100, 'trafo') ...
                             };

    testcases_online_trafo_exhaustive = { ...
      nfsftTestcaseDelegateOnline(128, 1000, 'trafo') ...
      nfsftTestcaseDelegateOnline(256, 1000, 'trafo') ...
      nfsftTestcaseDelegateOnline(512, 10, 'trafo') ...
      nfsftTestcaseDelegateOnline(512, 100, 'trafo') ...
      nfsftTestcaseDelegateOnline(1024, 10, 'trafo') ...
      nfsftTestcaseDelegateOnline(2048, 10, 'trafo') ...
    };

    testcases_online_adjoint = { ...
      nfsftTestcaseDelegateOnline(32, 10, 'adjoint') ...
      nfsftTestcaseDelegateOnline(32, 100, 'adjoint') ...
      nfsftTestcaseDelegateOnline(32, 1000, 'adjoint') ...
      nfsftTestcaseDelegateOnline(64, 10, 'adjoint') ...
      nfsftTestcaseDelegateOnline(64, 100, 'adjoint') ...
      nfsftTestcaseDelegateOnline(128, 10, 'adjoint') ...
      nfsftTestcaseDelegateOnline(128, 100, 'adjoint') ...
      nfsftTestcaseDelegateOnline(256, 10, 'adjoint') ...
      nfsftTestcaseDelegateOnline(256, 100, 'adjoint') ...
                       };

    testcases_online_adjoint_exhaustive = { ...
      nfsftTestcaseDelegateOnline(128, 1000, 'adjoint') ...
      nfsftTestcaseDelegateOnline(256, 1000, 'adjoint') ...
      nfsftTestcaseDelegateOnline(512, 10, 'adjoint') ...
      nfsftTestcaseDelegateOnline(512, 100, 'adjoint') ...
      nfsftTestcaseDelegateOnline(1024, 10, 'adjoint') ...
      nfsftTestcaseDelegateOnline(2048, 10, 'adjoint') ...
    };

    initializers_std = { ...
      nfsftTestcaseInitDelegate('init') ...
      nfsftTestcaseInitDelegate('init_advanced') ...
      nfsftTestcaseInitDelegate('init_guru',NFSFT_USE_DPT,FPT_NO_FAST_ALGORITHM,bitshift(1, 4),6) ...
      nfsftTestcaseInitDelegate('init_class') ...
      nfsftTestcaseInitDelegate('init_class', NFSFT_USE_DPT, FPT_NO_FAST_ALGORITHM) ...
    };

    testcases_online_trafo_equispaced = { ...
      nfsftTestcaseDelegateOnline(16, 0, 'trafo') ...
      nfsftTestcaseDelegateOnline(32, 0, 'trafo') ...
      nfsftTestcaseDelegateOnline(64, 0, 'trafo') ...
    };

    testcases_online_adjoint_equispaced = { ...
      nfsftTestcaseDelegateOnline(16, 0, 'adjoint') ...
      nfsftTestcaseDelegateOnline(32, 0, 'adjoint') ...
      nfsftTestcaseDelegateOnline(64, 0, 'adjoint') ...
                             };
  
    initializers_equispaced = { ...
      nfsftTestcaseInitDelegate('init_class',NFSFT_EQUISPACED) ...
    };

    testcases_gl_quadrature_online = { ...
      nfsftTestcaseDelegateGLQuadratureOnline(4) ...
      nfsftTestcaseDelegateGLQuadratureOnline(8) ...
      nfsftTestcaseDelegateGLQuadratureOnline(16) ...
      nfsftTestcaseDelegateGLQuadratureOnline(32) ...
      nfsftTestcaseDelegateGLQuadratureOnline(64) ...
      nfsftTestcaseDelegateGLQuadratureOnline(128) ...
      nfsftTestcaseDelegateGLQuadratureOnline(256) ...
      nfsftTestcaseDelegateGLQuadratureOnline(512) ...
    };

    testcases_gl_quadrature_online_large = { ...
      nfsftTestcaseDelegateGLQuadratureOnline(1024) ...
    };

    testcases_gl_quadrature_online_exhaustive = { ...
      nfsftTestcaseDelegateGLQuadratureOnline(2048) ...
    };
  
    initializers_quadrature = { ...
      nfsftTestcaseInitDelegate('init_advanced',NFSFT_NORMALIZED) ...
      nfsftTestcaseInitDelegate('init_advanced',bitor(NFSFT_NORMALIZED,NFSFT_USE_DPT)) ...
      nfsftTestcaseInitDelegate('init_class',NFSFT_NORMALIZED) ...
      nfsftTestcaseInitDelegate('init_class',bitor(NFSFT_NORMALIZED,NFSFT_USE_DPT)) ...
    };

    initializers_quadrature_large = { ...
      nfsftTestcaseInitDelegate('init_class',NFSFT_NORMALIZED) ...
    };
  end
  
  methods
    function ok = check_single(~, testcase, init_delegate, check_delegate, trafo_delegate)
      ok = 0;
      
      testcase = testcase.setup;
      
      if isempty(ver('octave'))
        bound = trafo_delegate.acc(init_delegate.nfft_cutoff_m);
        X = testcase.x;
        N = testcase.N;
        switch trafo_delegate.name
          case 'trafo'
            fh_testcase = testcase.f_hat;
            f_mat = legendre(0,cos(X(2,:)),'norm').' * fh_testcase(0,0) * sqrt(2);
            for n=1:N
              Lp = legendre(n,cos(X(2,:)),'norm') * sqrt(2/(2*n+1));
              f_mat=f_mat + (flipud(Lp).* exp((-n:0)'*1i*X(1,:))).'* fh_testcase(n,-n:0)...
                + (Lp(2:end,:).* exp((1:n)'*1i*X(1,:))).'* fh_testcase(n,1:n);
            end
            f_testcase = testcase.f;
            err = max(abs(f_testcase-f_mat)) / sum(sum(abs(double(fh_testcase))));
            if err > bound
              ok = 0;
              fprintf(' -> %-4s (incorrect reference val)\n','FAILED');
              return;
            end
          case 'adjoint'
            fh_testcase = testcase.f_hat;
            f_testcase = testcase.f;
            fh2_mat = zeros(size(testcase.f_hat));
            fh2_mat(1)= legendre(0,cos(X(2,:)),'norm') * sqrt(2) * f_testcase;
            for n=1:N
              Lp = legendre(n,cos(X(2,:)),'norm') * sqrt(2/(2*n+1));
              fh2_mat(n^2+1:n^2+n+1)= (flipud(Lp).* exp(-(-n:0)'*1i*X(1,:)))* f_testcase;
              fh2_mat(n^2+n+2:n^2+2*n+1)= (Lp(2:end,:).* exp(-(1:n)'*1i*X(1,:)))* f_testcase;
            end
            fh2_mat = f_hat(fh2_mat);
            err = max(max(abs(double(fh_testcase-fh2_mat)))) / sum(abs(f_testcase));
            if err > bound
              ok = 0;
              fprintf(' -> %-4s (incorrect reference val)\n','FAILED');
              return;
            end
        end
      end

      init_txt = init_delegate.txt;
      fprintf(', %-18s', init_txt);
      [init_delegate, plan] = init_delegate.init(testcase.N, testcase.M);

      if ~isempty(init_delegate.nfft_cutoff_m)
        fprintf(', m = %2d', init_delegate.nfft_cutoff_m);
      else
        fprintf(', m std.');
      end
      fprintf(', %-14s', trafo_delegate.name);

      % Equispaced nodes are automatically set in nfsft_init
      if ~bitand(init_delegate.flags, NFSFT_EQUISPACED)
        if isnumeric(plan)
          nfsft_set_x(plan, testcase.x);
          nfsft_precompute_x(plan);
        else
          plan.x = testcase.x;
        end
      end

      run_test = 1;
      if (trafo_delegate.isCost)
        cost = trafo_delegate.cost(plan);
        if (cost > 0.1)
          fprintf(' -> %-4s (cost too high)\n','OK');
          ok = 1;
          run_test = 0;
        end
      end
      
      if run_test
        plan = check_delegate.prepare(plan, testcase.f, testcase.f_hat);
        
        plan = trafo_delegate.trafo(plan);
        
        err = check_delegate.compare(plan, testcase.f, testcase.f_hat);
        bound = trafo_delegate.acc(init_delegate.nfft_cutoff_m);
        ok = err < bound;
        if ~ok; okstr = 'FAIL'; else, okstr = 'OK'; end
        fprintf(' -> %-4s %.3e (%.3e)\n', okstr, err, bound);
      end
      
      testcase.destroy;
      if isnumeric(plan)
        nfsft_finalize(plan);
      else
        delete(plan);
      end
      nfsft_forget;
    end

    function result = check_many(h, testcases, initializers, check_delegate, trafos)
      result = 1;
      for k=1:length(trafos)
        for i=1:length(testcases)
          for j=1:length(initializers)
            r = h.check_single(testcases{i}, initializers{j}, check_delegate, trafos{k});
            result = min(result, r);
          end
        end
      end
    end

    function result = nfsft_check_online(h)
      result = h.check_many(h.testcases_online_trafo, h.initializers_std, nfsftTestcaseCheckDelegate('trafo'), {nfsftTestcaseTrafoDelegate('trafo')});
      if h.perform_exhaustive_tests_flag
        result2 = h.check_many(h.testcases_online_trafo_exhaustive, h.initializers_std, nfsftTestcaseCheckDelegate('trafo'), {nfsftTestcaseTrafoDelegate('trafo')});
        result = min(result,result2);
      end
    end
    
    function result = nfsft_check_adjoint_online(h)
      result = h.check_many(h.testcases_online_adjoint, h.initializers_std, nfsftTestcaseCheckDelegate('adjoint'), {nfsftTestcaseTrafoDelegate('adjoint')});
      if h.perform_exhaustive_tests_flag
        result2 = h.check_many(h.testcases_online_adjoint_exhaustive, h.initializers_std, nfsftTestcaseCheckDelegate('adjoint'), {nfsftTestcaseTrafoDelegate('adjoint')});
        result = min(result,result2);
      end
    end

    function result = nfsft_check_trafo_equispaced_online(h)
      result = h.check_many(h.testcases_online_trafo_equispaced, h.initializers_equispaced, nfsftTestcaseCheckDelegate('trafo'), {nfsftTestcaseTrafoDelegate('trafo')});
    end

    function result = nfsft_check_adjoint_equispaced_online(h)
      result = h.check_many(h.testcases_online_adjoint_equispaced, h.initializers_equispaced, nfsftTestcaseCheckDelegate('adjoint'), {nfsftTestcaseTrafoDelegate('adjoint')});
    end

    function ok = check_single_quadrature(~, testcase, init_delegate, trafo_delegate)
      ok = 0;
      
      init_txt = init_delegate.txt;
      fprintf('... %-18s', init_txt);
      [init_delegate, plan] = init_delegate.init(testcase.N, testcase.M);

      if ~isempty(init_delegate.nfft_cutoff_m)
        fprintf(', m = %2d', init_delegate.nfft_cutoff_m);
      else
        fprintf(', m std.');
      end
      fprintf(', %-14s', trafo_delegate.name);

      if isnumeric(plan)
        nfsft_set_x(plan, testcase.x);
        nfsft_precompute_x(plan);
      else
        plan.x = testcase.x;
      end

      run_test = 1;
      if (trafo_delegate.isCost)
        cost = trafo_delegate.cost(plan);
        if (cost > 0.1)
          fprintf(' -> %-4s (cost too high)\n','OK');
          ok = 1;
          run_test = 0;
        end
      end
      
      if run_test
        fh = testcase.f_hat;

        f_mul_W = testcase.f;
        if isnumeric(plan)
          nfsft_set_f(plan,f_mul_W);
        else
          plan.f = f_mul_W;
        end
          
        plan = trafo_delegate.trafo(plan);

        if isnumeric(plan)
          fh2 = nfsft_get_f_hat(plan);
        else
          fh2 = plan.fhat;
        end

        err = norm(fh-fh2) / norm(fh);
        bound = 2^22* eps;
        ok = err < bound;
        if ~ok; okstr = 'FAIL'; else, okstr = 'OK'; end
        fprintf(' -> %-4s %.3e (%.3e)\n', okstr, err, bound);
      end
      
      if isnumeric(plan)
        nfsft_finalize(plan);
      else
        delete(plan);
      end
    end
    
    function result = check_many_quadrature(h, testcases, initializers, trafos)
      result = 1;
      for i=1:length(testcases)
        testcase = testcases{i}.setup;
        bound = 2^22* eps;
        err = max(abs(testcase.f-testcase.f2)) / sum(sum(abs(double(testcase.f_hat))));
        ok = err < bound;
        if ~ok
          fprintf('Testcase setup FAILED\n');
          result = 0;
          continue;
        end

        fprintf('\n');

        for j=1:length(initializers)
          for k=1:length(trafos)
            r = h.check_single_quadrature(testcase, initializers{j}, trafos{k});
            result = min(result, r);
          end
          nfsft_forget;
        end
        testcase.destroy;
      end
    end

    function result = nfsft_check_quadrature_online(h)
      result = h.check_many_quadrature(h.testcases_gl_quadrature_online, h.initializers_quadrature, {nfsftTestcaseTrafoDelegate('adjoint')});

      % Don't perform large test in 32bit Octave
      archstring = computer;
      if isempty(strfind(archstring,'i686'))
        result2 = h.check_many_quadrature(h.testcases_gl_quadrature_online_large, h.initializers_quadrature_large, {nfsftTestcaseTrafoDelegate('adjoint')});
        result = min(result,result2);
      end

      if h.perform_exhaustive_tests_flag
        result2 = h.check_many_quadrature(h.testcases_gl_quadrature_online_exhaustive, h.initializers_quadrature, {nfsftTestcaseTrafoDelegate('adjoint')});
        result = min(result,result2);
      end
    end

  end
end