function gridpts = generate_Smolyak_grid(d,w,varargin)
% generates a Smolyak grid in d-dimensions, with level w
% for [-1,1]^d
% default: Clenshaw-Curtis grid
% grid types: 1 - Clenshaw-Curtis
%             2 - Gaussian


if isempty(varargin)
    grid_type = 1;
else
    grid_type = varargin{1};
end
if grid_type == 1
    ini_idx = w+d;
else
    ini_idx = w+1;
end

gridpts = [];

fin_idx = w+d;
for idx=ini_idx:fin_idx
    i_idxs = form_integer_partition(idx,d); %indices of the i in the partitioning
    nPartitions = size(i_idxs,1);
    collect_abs = cell(d,1); %collection of abscissa in each dimension

    % for each partition, obtain the corresponding quadrature points
    for i=1:nPartitions
        curr_partition = i_idxs(i,:); 
        curr_m = 2 .^ (curr_partition-1) + 1; % this is the clenshaw-curtis nodes
        curr_m (curr_partition == 1) = 1;

        %access the abscissa points in each dimension, can either be based
        %on Clenshaw-Curtis quad or Gaussian quad
        for j=1:d
            dim_pts = curr_m(j); %the number of points in the current dimension (for the current partition)
            collect_abs{j} = zeros(1,dim_pts);
            if dim_pts == 1
                abs_points = 0;
            else
                if grid_type == 1
                    abs_points = -cos(pi*((1:dim_pts)-1)/(dim_pts-1)); %  Clenshaw-curtis nodes
                else
                    abs_points = lgwt(dim_pts,-1,1)'; %  Legendre-Gauss nodes
                end
            end
            collect_abs{j} = abs_points;
        end
        [abs_nowCollect{1:d}] = ndgrid(collect_abs{1:d}); % build tensor product
        nGridPts = prod(curr_m);
        s_now = zeros(d,nGridPts);
        for j=1:d
            s_now(j,:) = abs_nowCollect{j}(:)';
        end
        gridpts = [gridpts,s_now]; % stacking the grids
    end
end


%plots the Smolyak grid points in dimensions 2 or 3
figure();
if d==2
    grid_x = gridpts(1,:);
    grid_y = gridpts(2,:);
    plot(grid_x(:), grid_y(:), 'r.', 'MarkerSize', 6);
    axis equal;
    xlim([-1 1]); ylim([-1 1]);
    box on;
%     plot(grid_x,grid_y,"*");
elseif d==3
    grid_x = gridpts(1,:);
    grid_y = gridpts(2,:);
    grid_z = gridpts(3,:);
    plot3(grid_x,grid_y,grid_z,"*");
end
