def tree(X, n, symbolic=False):
    """Evaluates the entire tree of orthogonal polynomials for the n-cube

    The computation is organized such that tree returns a list of arrays, L={0,
    ..., dim}, where each level corresponds to the polynomial degree L.
    Further, each level is organized like a discrete (dim-1)-dimensional
    simplex. Let's demonstrate this for 3D:

    L = 1:
         (0, 0, 0)

    L = 2:
         (1, 0, 0)
         (0, 1, 0) (0, 0, 1)

    L = 3:
         (2, 0, 0)
         (1, 1, 0) (1, 0, 1)
         (0, 2, 0) (0, 1, 1) (0, 0, 2)

    The main insight here that makes computation for n dimensions easy is that
    the next level is composed by:

       * Taking the whole previous level and adding +1 to the first entry.
       * Taking the last row of the previous level and adding +1 to the second
         entry.
       * Taking the last entry of the last row of the previous and adding +1 to
         the third entry.

    In the same manner this can be repeated for `dim` dimensions.
    """
    p0, a, b, c = legendre(n + 1, "normal", symbolic=symbolic)

    dim = X.shape[0]

    p0n = p0 ** dim
    out = []

    level = numpy.array([numpy.ones(X.shape[1:], dtype=int) * p0n])
    out.append(level)

    # TODO use a simpler binom implementation
    for L in range(n):
        level = []
        for i in range(dim - 1):
            m1 = int(scipy.special.binom(L + dim - i - 1, dim - i - 1))
            if L > 0:
                m2 = int(scipy.special.binom(L + dim - i - 2, dim - i - 1))
            r = 0
            for k in range(L + 1):
                m = int(scipy.special.binom(k + dim - i - 2, dim - i - 2))
                val = out[L][-m1:][r : r + m] * (a[L - k] * X[i] - b[L - k])
                if L - k > 0:
                    val -= out[L - 1][-m2:][r : r + m] * c[L - k]
                r += m
                level.append(val)

        # treat the last one separately
        val = out[L][-1] * (a[L] * X[-1] - b[L])
        if L > 0:
            val -= out[L - 1][-1] * c[L]
        level.append([val])

        out.append(numpy.concatenate(level))

    return out
def tree ( X , n , symbolic = False ) : p0 , a , b , c = legendre ( n + 1 , "normal" , symbolic = symbolic ) dim = X . shape [ 0 ] p0n = p0 ** dim out = [ ] level = numpy . array ( [ numpy . ones ( X . shape [ 1 : ] , dtype = int ) * p0n ] ) out . append ( level ) for L in range ( n ) : level = [ ] for i in range ( dim - 1 ) : m1 = int ( scipy . special . binom ( L + dim - i - 1 , dim - i - 1 ) ) if L > 0 : m2 = int ( scipy . special . binom ( L + dim - i - 2 , dim - i - 1 ) ) r = 0 for k in range ( L + 1 ) : m = int ( scipy . special . binom ( k + dim - i - 2 , dim - i - 2 ) ) val = out [ L ] [ - m1 : ] [ r : r + m ] * ( a [ L - k ] * X [ i ] - b [ L - k ] ) if L - k > 0 : val -= out [ L - 1 ] [ - m2 : ] [ r : r + m ] * c [ L - k ] r += m level . append ( val ) val = out [ L ] [ - 1 ] * ( a [ L ] * X [ - 1 ] - b [ L ] ) if L > 0 : val -= out [ L - 1 ] [ - 1 ] * c [ L ] level . append ( [ val ] ) out . append ( numpy . concatenate ( level ) ) return out
Evaluates the entire tree of orthogonal polynomials for the n-cube

    The computation is organized such that tree returns a list of arrays, L={0,
    ..., dim}, where each level corresponds to the polynomial degree L.
    Further, each level is organized like a discrete (dim-1)-dimensional
    simplex. Let's demonstrate this for 3D:

    L = 1:
         (0, 0, 0)

    L = 2:
         (1, 0, 0)
         (0, 1, 0) (0, 0, 1)

    L = 3:
         (2, 0, 0)
         (1, 1, 0) (1, 0, 1)
         (0, 2, 0) (0, 1, 1) (0, 0, 2)

    The main insight here that makes computation for n dimensions easy is that
    the next level is composed by:

       * Taking the whole previous level and adding +1 to the first entry.
       * Taking the last row of the previous level and adding +1 to the second
         entry.
       * Taking the last entry of the last row of the previous and adding +1 to
         the third entry.

    In the same manner this can be repeated for `dim` dimensions.