Total degree basis

[heptaflar]

Example:

import itertools

D = 10  # Dimensions
N = 3   # Maximum Total Degree

# Generate all combinations of indices from 0 to N for 10 dimensions
all_combinations = itertools.product(range(N + 1), repeat=D)

# Filter: ONLY keep combinations where the SUM of indices is <= N
total_degree_indices = [idx for idx in all_combinations if sum(idx) <= N]

print(f"Total coefficients kept: {len(total_degree_indices)}")
# For N=3, D=10, this prints 286 (instead of 4^10 = 1,048,576!)

[heptaflar]

import jax
import jax.numpy as jnp
from typing import Tuple

# Assuming your b_basis function looks something like this:
# def b_basis(deg, x_val): ... 

@jax.jit
def b_poly_point_total_degree(c: jnp.ndarray, x: jnp.ndarray, degrees: Tuple[int, ...], indices: jnp.ndarray) -> jnp.ndarray:
    """
    Evaluates an N-variate Total Degree Bernstein polynomial on a single point.

    :param c: Flat 1D array of Bernstein coefficients, shape (M,)
    :param x: The point x, shape (N,)
    :param degrees: A tuple of max degrees for each dimension, e.g., (3, 3, ..., 3)
    :param indices: Precomputed 2D array of valid Total Degree indices, shape (M, N)
                    where sum(indices[i]) <= max_total_degree
    :returns: The polynomial value B(x) as a scalar/0D array.
    """
    N = x.shape[0]
    
    # 1. Compute 1D basis values for ALL possible degrees across all dimensions
    # To keep JAX happy with static shapes, we find the absolute max degree
    max_deg = max(degrees)
    
    # Pre-evaluate the 1D basis functions for each dimension
    # We build a matrix 'basis_matrix' of shape (N, max_deg + 1)
    # where basis_matrix[dim, k] is the k-th Bernstein basis value for x[dim]
    
    # Note: If your b_basis function handles vectors, you can vectorize this.
    # For simplicity and JAX compliance, we can stack them:
    basis_list = [b_basis(degrees[dim], x[dim]) for dim in range(N)]
    
    # Pad lists with zeros if dimensions have different max degrees, then stack
    # If all dimensions have the same degree, it's just: jnp.stack(basis_list)
    basis_matrix = jnp.stack([jnp.pad(b, (0, max_deg + 1 - len(b))) for b in basis_list])

    # 2. For each active coefficient, look up its corresponding 1D basis values
    # indices shape: (M, N). We map over the N dimensions.
    # row_indices creates a grid to pluck out the exact basis value for each dimension
    dim_selectors = jnp.arange(N)
    
    # Evaluated terms shape: (M, N) -> holds the 1D basis value for each variable in each term
    evaluated_terms = basis_matrix[dim_selectors, indices]
    
    # 3. Multiply across the dimensions to get the full multi-variate basis value for each term
    # product_terms shape: (M,)
    product_terms = jnp.prod(evaluated_terms, axis=1)

    # 4. Final dot product: sum(c_i * basis_i)
    return jnp.dot(c, product_terms)