function [res,err] = doubly_adaptive(f,x,varargin)

%%%%%
%
% INPUT: 
%   f:      Function to integrate
%   x:      Region 
% 
%   Named Arguments: 
%   tol:    Tolerance
%   F:      Hint parameter
%   maxdiv: Max number of domain subdivisions
%   maxord: Maximum order for the formulas
%   err_method: Error estimation method
%           1: Higher order Gauss
%           2: Antigauss
%           3: Gauss-Kronrod
%   domain:
%           1: Torus
%           2: Sphere
%   coord: Coordinate system
%           0: Cartesian
%           1: Polar (spherical or toroidal)
%
% OUTPUT:
%   res:    Computer integral
%   err:    Error estimate
%
%%%%%

%# define defaults at the beginning of the code so that you do not need to
%# scroll way down in case you want to change something or if the help is
%# incomplete
options = struct('tol',10^(-15),'f',0.1,'maxdiv',5,'maxord',20,'err_method',1,'domain',1,'coord',0);

%# read the acceptable names
optionNames = fieldnames(options);

%# count arguments
nArgs = length(varargin);
if round(nArgs/2)~=nArgs/2
   error('EXAMPLE needs propertyName/propertyValue pairs')
end

for pair = reshape(varargin,2,[]) %# pair is {propName;propValue}
   inpName = lower(pair{1}); %# make case insensitive

   if any(strmatch(inpName,optionNames))
      %# overwrite options. If you want you can test for the right class here
      %# Also, if you find out that there is an option you keep getting wrong,
      %# you can use "if strcmp(inpName,'problemOption'),testMore,end"-statements
      options.(inpName) = pair{2};
   else
      error('%s is not a recognized parameter name',inpName)
   end
end

err_method = options.('err_method');
maxdiv = options.('maxdiv');
maxord = options.('maxord');
tol = options.('tol');
F = options.('f');
domain = options.('domain');
coord = options.('coord');

%tic

cubature = cellfun(@str2func,{'cubature_2_torus' 'cubature_2_subsphere_fast'},'UniformOutput',false);
counter = 0;

% Init some variables;
queue = [];
done = [];

% Some matrix juggling because otherwise Matlab uses only first return value
x = [x(1,:),x(2,:)];
[res,err,hint,ord,div,reg,evals] = init_region(f,x,0,err_method,domain,coord);
counter = counter + evals;
% This will be used when we split the region
queue = [queue; res,err,hint,ord,div,reg];

res = queue(end,1);
err = queue(end,2);

% Stop when error goes below threshold in input or we've worked on all regions
while (size(queue,1) > 0 && ...
    err > tol               ...
    )
    % Work on the region with highest error
    [m,index] = max(queue(:,2));
    rect = queue(index,:);
    
    % If the hint is below F, try order adaptive algorithm
    if rect(3) < F
        %disp('Going to higher order');

        % If we've reached maximum order, move this region from
        % the original queue to the "done" one
        if rect(4) >= maxord
            queue = [queue(1:index-1,:); queue(index+1:end,:)];
            done = [done; rect];

        % Else, update the info on this region in the queue
        % using higher order formula
        else
            [tw,xw] = cubature{domain}([rect(6:7);rect(8:end)] ...
            ,rect(4)+1,0,1);

            %fx = feval(f,xw(:,1),xw(:,2),xw(:,3));
            if coord > 0
                fx=feval(f,tw(:,1),tw(:,2));
                w=tw(:,3);
            else
                fx=feval(f,xw(:,1),xw(:,2),xw(:,3));
                w=xw(:,4);
            end
            counter = counter + 1;
            %w = xw(:,4);
            p_res = w'*fx;

            h_res = rect(1);
            if err_method > 1
                % Higher order rule for error estimate
                [tw,xw] = cubature{domain}([rect(6:7);rect(8:end)] ...
                ,rect(4)+1,0,err_method);

                fx = feval(f,xw(:,1),xw(:,2),xw(:,3));
                w = xw(:,4);
                h_res = w'*fx;
                % End of higher order calculations
            end
            %printf('Valutazioni della funzione f: %d\n',counter);
            p_err = abs(h_res - p_res);
            if err_method == 2
                p_err = p_err/2;
            end
            rect(4) = rect(4) * 2;
            rect(3) = p_err/rect(2);
            rect(2) = p_err;
            rect(1) = p_res;
            queue = [queue(1:index-1,:); rect; queue(index+1:end,:)];
        end
    % Else, try splitting the domain
    else

        % If we've already split this region enough times,
        % move it to the "done" queue
        if rect(5) >= maxdiv
            queue = [queue(1:index-1,:); queue(index+1:end,:)];
            done = [done; rect];

        % Else, split it and add the new regions to the queue
        else
            xs(:,1) = rect(6:end);
            xs(:,2) = rect(6:end);
            xs(:,3) = rect(6:end);
            xs(:,4) = rect(6:end);
            hdist = (xs(2,1)-xs(1,1))/2;
            vdist = (xs(4,1)-xs(3,1))/2;

            xs(2,1) = xs(2,1) - hdist;
            xs(4,1) = xs(4,1) - vdist;
            xs(1,2) = xs(1,2) + hdist;
            xs(4,2) = xs(4,2) - vdist;
            xs(2,3) = xs(2,3) - hdist;
            xs(3,3) = xs(3,3) + vdist;
            xs(1,4) = xs(1,4) + hdist;
            xs(3,4) = xs(3,4) + vdist;
            part_queue = [queue(1:index-1,:)];
            for j=1:4
                [res,err,hint,ord,div,reg,evals] = init_region(f,xs(:,j)',rect(5),err_method,domain,coord);
                counter = counter + evals;
                part_queue = [part_queue; res, err, hint, ord, div, reg];
            end
            part_queue = [part_queue; queue(index+1:end,:)];
            queue = part_queue;
        end
    end
    % Update new result and error
    % so we can check the "while" condition
    res = 0;
    err = 0;
    if size(queue,2) >= 2
        res  = res + sum(queue(:,1));
        err  = err + sum(queue(:,2));
    end
    if size(done,2) >= 2
        res  = res + sum(done(:,1));
        err  = err + sum(done(:,2));
    end
end

printf('FINAL_NONOPT: Valutazioni della funzione f: %d\n',counter);
%toc
end

function [res,err,hint,ord,div,x,p_evals] = init_region(f,x,div,err_method,domain,coord)

%%%%%
% Generate data structure for new regions
% This data will be the initial queue element for the region
%
% Input: integrand, region, and times it was already split
% Output: integral, error estimate, hint, order, 
% times the region was split and the region itself
%%%%%
%disp('Init-ing new domain');
init = zeros(4,2);
cubature = cellfun(@str2func,{'cubature_2_torus' 'cubature_2_subsphere_fast'},'UniformOutput',false);
p_evals = 0;

for i=2:4
    [tw,xw] = cubature{domain}([x(1:2);x(3:end)],i,0,1);

    if coord > 0
        fx=feval(f,tw(:,1),tw(:,2));
        w=tw(:,3);
    else
        fx=feval(f,xw(:,1),xw(:,2),xw(:,3));
        w=xw(:,4);
    end
    p_evals = p_evals + 1;
    %w = xw(:,4);
    init(i,1) = w'*fx;
    h_res = init(i-1,1);
    if err_method > 1
        % Higher order rule for error estimate
        [tw,xw] = cubature{domain}([x(1:2);x(3:end)] ...
        ,i,0,err_method);

        fx = feval(f,xw(:,1),xw(:,2),xw(:,3));
        p_evals = p_evals + 1;
        w = xw(:,4);
        h_res = w'*fx;
        % End of higher order calculations
    end
    init(i,2) = abs(init(i,1) - h_res);
    if err_method == 2
        init(i,2) = init(i,2)/2;
    end
end

res = init(4,1);
err = init(4,2);
ord = i;
div = div+1;
hint = init(4,2)/init(3,2);
end
