# Smolyak_Elem_Isotrop.m is a routine that constructs the subindices of the
# Smolyak elements (grid points and basis functions) for the isotropic case;
# see "Smolyak method for solving dynamic economic models: Lagrange interpo-
# lation, anisotropic grid and adaptive domain" by Kenneth L. Judd, Lilia Maliar,
# Serguei Maliar and Rafael Valero, (2014), Journal of Economic Dynamics and
# Control 44, 92�123 (henceforth, JMMV (2014)), Section 2.2.3
#
# This version: Novenber 5, 2014. First version: December 17, 2012.
# -------------------------------------------------------------------------
# Input:   "d" is the number of dimensions (the number of state  variables)
#          "mu" is the level of approximation
#
# Output:  "Smolyak_elem_iso" is the vector of subindices of unidimensional
#           elements (Smolyak grid points or polynomial basis function) that
#           constitute a multidimensional element (Smolyak grid point or
#           polynomial basis function) for the isotropic case
#
# -------------------------------------------------------------------------
# Copyright � 2014 by Lilia Maliar, Serguei Maliar and Rafael Valero. All
# rights reserved. The code may be used, modified and redistributed under
# the terms provided in the file "License_Agreement.txt".
# -------------------------------------------------------------------------

 function Smolyak_Elem_Isotrop(d,mu)

# 1. Identify the indices of disjoint sets A's, i1,...,id, that satisfy the
# Smolyak rule, d<=i1+i2+...+id<=|i|; see equation (1) in JMMV(2014)
# -------------------------------------------------------------------------

    Smol_rule = [];     # This will be a matrix of subindices i1,...,id of
                        # disjoint sets A_i's such that their sum across all
                        # dimensions, |i|=i1+i2+...+id,  satisfies  the
                        # Smolyak rule, d<=|i|<=d+mu; initially, the matrix
                        # is empty; at the intermediate steps j=0,...,mu,
                        # the subindices satisfy d<=|i|=d+j

    incr_Smol_rule = [];# The matrix of subindices i1,...,id of disjoint
                        # sets A_i's, such that their sum across all dimensions,
                        # |i|=i1+i2+...+id, is exactly equal to d+j, i.e.,
                        # |i|=d+j; initially, the matrix is empty; when j
                        # increases, this matrix is concatinated to matrix
                        # "Smol_rule" obtained for the previous value of j,
                        # i.e., j-1
   for j = 0:mu

     prev_incr = incr_Smol_rule;
                       # Ffor the previous value of j, call "incr_Smol_rule"
                       # as "prev_incr"

     # Identify new subindices of unidimensional sets A's i1,i2,...,id that
     # jointly satisfy the Smolyak rule as the sum of subindices increases
     # from j-1 to j; see JMMV (2014), Section 2.2
     #---------------------------------------------------------------------
     if (j == 0)
       incr_Smol_rule = ones(1,d); # The matrix of subindices for j=0 is a
                                   # 1-by-d vector of ones, (1,...,1)

     else
       m = size(prev_incr,1);
                                   # Compute the number of columns in the
                                   # previous matrix of subindices that
                                   # satisfy the Smolyak rule "prev_incr"
       incr_Smol_rule = [];        # Initially, the matrix of subindices is
                                   # empty
       #aux = zeros(m,d);           # Allocate memory to an initial auxiliary
                                   # matrix that will be added to
                                   # "prev_incr"
        for id = 1:d
            aux_new = zeros(m,d);
            #aux_new = aux;         # New auxiliary matrix is equal to the old
                                   # one
            aux_new[:,id] .= 1;     # For a dimension i, set elements of this
                                   # new auxiliary matrix to 1 instead of 0

            augmented = prev_incr + aux_new;
                                   # Increase the subinices of
                                   # "prevoius_incr" by 1 in dimension id

            if size(incr_Smol_rule,1)==0
                incr_Smol_rule=augmented
            else
            incr_Smol_rule = vcat(incr_Smol_rule, augmented);
                                   # Concatenate "incr_Smol_rule" and
                                   # "augmented": the first row of "augmented"
                                   # goes after the last row of
                                   # "incr_Smol_rule"
            end
        end


     end
         incr_Smol_rule = sortslices(unique(incr_Smol_rule,dims=1),dims=1);
                                   # Eliminate the repeated indices in the
                                   # matrix "incr_Smol_rule"

     # Concatenate the matrix of newly constructed indices to the previously
     # obtained matrix of indices satisfying the Smolyak rule
     #---------------------------------------------------------------------
     if size(Smol_rule,1)==0
         Smol_rule=incr_Smol_rule;
     else
     Smol_rule = vcat(Smol_rule, incr_Smol_rule); # E.g., for mu=1 and d=2,
                                                   # Smol_rule=[1 1; 1 2; 2 1]
     end
   end

   n_comb = size(Smol_rule,1);
              # The total number of combinations of indices, i1,...,id,
              # of the unidimensional disjoint sets A's that satisfy the
              # Smolyak rule; e.g., for mu=1 and d=2, n_comb=3 such as (1)
              # i1=1 and i2=1, (2) i1=1 and i2=2, (3) i1=2 and i2=1


# 2. Construct the multidimensional indices of elements as a Cartesian product
# of unidimensional indices
# -------------------------------------------------------------------------

Smolyak_elem_iso = [];
                # The matrix of multidimensional indices of elements (points)
                #  belonging to the disjoint subsets A's that satisfy the
                # Smolyak rule (every single point is indexed); initially,
                # the matrix is empty

for i = 1:n_comb            # For each combination of subindices i1,...,id
                            # that satisfies the Smolyak rule


    incr_indices = [];      # This is a matrix of multidimensional indices
                            # of unidimensional grid points that will be
                            # added to the final multidimensional matrix of
                            # indices "Smolyak_elem_iso", as n_comb increases;
                            # initially, the matrix is empty

    one_comb = Smol_rule[i,:];
                            # Consider the i-th row of "Smol_rule" (i.e.,
                            # one particular combination of the subindices
                            # of the disjoint sets satisfying the Smolyak
                            # rule); e.g., for mu=1 and d=2, Smol_rule=
                            # [1 1; 1 2; 2 1] and one_comb for i=1 is [1 1];
                            # 1-by-d
    for jd = 1:d            # For each dimension jd, ...

      prev_indices = incr_indices;

          # Compute the indices of elements (points) of the unidimensional
          # set
          #----------------------------------------------------------------
            if one_comb[jd] == 1
                           # Take a jd-th element of the row vector "one_comb"
                           # that corresponds to dimension jd (this is the
                           # subindex of the unidimensional disjoint set A_i
                           # from which this element comes from; if an element
                           # (point) is from the disjoint set
                           # A_1,...
               indices_elem_jd = 1;
                           # A_1 contains one element; this element is indexed "1"

            elseif one_comb[jd] == 2
                           # If an element (point) is from the disjoint set
                           # A_2,...
               indices_elem_jd = [2;2^(one_comb[jd]-1)+1];
                           # A_2 contains two elements; these elements are
                           # indexed "2" and "3", so that indices_elem_jd=[2;3]
            else
               indices_elem_jd = collect(2^(one_comb[jd]-2)+2 : 2^(one_comb[jd]-1)+1);
                           # The subsequent disjoint sets contain the elements
                           # from m(one_comb(jd)-1)+1=2^(one_comb(jd)-2)+2 to
                           # m(one_comb(jd))=2^(one_comb(jd)-1)+1; e.g., the
                           # disjoint set A_4 contains the elements indexed
                           # 6, 7, 8 and 9, so that indices_elem_jd=[6,7,8,9]
            end

           # Create a Cartesian product of two sets, "prev_indices" and
           # "indices_elem_jd"
           #-----------------------------------------------------------
           a = prev_indices;   # Call one set "a"
           b = indices_elem_jd;# Call the other set "b"
           if (isempty(b))     # If "b" is empty, ...
               z = a;          # The Cartesian product "z" is given by "a"
           elseif (isempty(a)) # If "a" is empty, ...
               z = b;          # The Cartesian product "z" is given by "b"
           else                # For the case of non-empty sets, ...

               z = [];         # Initially, the Cartesian product is empty
               a_columns = size(a,1);
               b_columns = size(b,1);
               for k = 1:b_columns
                  if isempty(size(b))
                      b_temp=b[k]
                  else
                      b_temp=b[k,:]
                  end
                  if isempty(z)
                      z = hcat(a,ones(a_columns,1)*b_temp);
                  else
                  z = vcat(z,hcat(a,ones(a_columns,1)*b_temp));
                  end
               end
           end

           incr_indices = z;
    end
    if isempty(Smolyak_elem_iso)
        Smolyak_elem_iso = incr_indices;
    else
        Smolyak_elem_iso = vcat(Smolyak_elem_iso,incr_indices);
                           # Construct the matrix of multidimensional indices
                           # of elements by concatenating "incr_indices" to
                           # the previous "Smolyak_elem_iso" obtained for the
                           # previous combination "n_comb"
    end
end

return Smolyak_elem_iso
end
