How to use the chaospy.poly function in chaospy

def __div__(self, other):
        """Python2 division."""
        return chaospy.poly.collection.arithmetics.mul(
            self, numpy.asfarray(other)**-1)
        return NotImplemented

def _stieltjes_approx(dist, order, accuracy, normed, **kws):
    """Stieltjes' method with approximative recurrence coefficients."""
    kws["rule"] = kws.get("rule", "C")
    assert kws["rule"].upper() != "G"

    absisas, weights = chaospy.quad.generate_quadrature(
        accuracy, dist.range(), **kws)
    weights = weights*dist.pdf(absisas)

    poly = chaospy.poly.variable(len(dist))
    orth = [poly*0, poly**0]

    inner = numpy.sum(absisas*weights, -1)
    norms = [numpy.ones(len(dist)), numpy.ones(len(dist))]
    coeff1 = []
    coeff2 = []

    for _ in range(order):

        coeff1.append(inner/norms[-1])
        coeff2.append(norms[-1]/norms[-2])
        orth.append((poly-coeff1[-1])*orth[-1] - orth[-2]*coeff2[-1])

        raw_nodes = orth[-1](*absisas)**2*weights
        inner = numpy.sum(absisas*raw_nodes, -1)
        norms.append(numpy.sum(raw_nodes, -1))

>>> orth, norms, A, B = cp.stieltjes(dist, 2, retall=True)
>>> print(cp.around(orth[2], 8))
[q0^2-q0+0.16666667]
>>> print(norms)
[[ 1.          0.08333333  0.00555556]]
    """

    assert not dist.dependent()

    try:
        dim = len(dist)
        K = np.arange(order+1).repeat(dim).reshape(order+1, dim).T
        A,B = dist.ttr(K)
        B[:,0] = 1.

        x = cp.poly.variable(dim)
        if normed:
            orth = [x**0*np.ones(dim), (x-A[:,0])/np.sqrt(B[:,1])]
            for n in range(1,order):
                orth.append((orth[-1]*(x-A[:,n]) - orth[-2]*np.sqrt(B[:,n]))/np.sqrt(B[:,n+1]))
            norms = np.ones(B.shape)
        else:
            orth = [x-x, x**0*np.ones(dim)]
            for n in range(order):
                orth.append(orth[-1]*(x-A[:,n]) - orth[-2]*B[:,n])
            orth = orth[1:]
            norms = np.cumprod(B, 1)

    except NotImplementedError:

        bnd = dist.range()
        kws["rule"] = kws.get("rule", "C")

def ensure_dtype(core, dtype, dtype_):
    """Ensure dtype is correct."""
    core = core.copy()
    if dtype is None:
        dtype = dtype_

    if dtype_ == dtype:
        return core, dtype

    for key, val in {
            int: chaospy.poly.typing.asint,
            float: chaospy.poly.typing.asfloat,
            np.float32: chaospy.poly.typing.asfloat,
            np.float64: chaospy.poly.typing.asfloat,
    }.items():

        if dtype == key:
            converter = val
            break
    else:
        raise ValueError("dtype not recognised (%s)" % str(dtype))

    for key, val in core.items():
        core[key] = converter(val)
    return core, dtype

>>> print(cp.outer(x, P))
        [q0, q0^2, q0^3]
        >>> print(cp.outer(P, P))
        [[1, q0, q0^2], [q0, q0^2, q0^3], [q0^2, q0^3, q0^4]]
    """
    if len(args) > 2:
        part1 = args[0]
        part2 = outer(*args[1:])

    elif len(args) == 2:
        part1, part2 = args

    else:
        return args[0]

    dtype = chaospy.poly.typing.dtyping(part1, part2)

    if dtype in (list, tuple, numpy.ndarray):

        part1 = numpy.array(part1)
        part2 = numpy.array(part2)
        shape = part1.shape +  part2.shape
        return numpy.outer(
            chaospy.poly.shaping.flatten(part1),
            chaospy.poly.shaping.flatten(part2),
        )

    if dtype == Poly:

        if isinstance(part1, Poly) and isinstance(part2, Poly):

            if (1,) in (part1.shape, part2.shape):

polynomial array, the output is the Jacobian matrix.

    Examples:
        >>> q0, q1 = chaospy.variable(2)
        >>> poly = chaospy.Poly([1, q0, q0*q1**2+1])
        >>> print(poly)
        [1, q0, q0q1^2+1]
        >>> print(differential(poly, q0))
        [0, 1, q1^2]
        >>> print(differential(poly, q1))
        [0, 0, 2q0q1]
    """
    poly = Poly(poly)
    diffvar = Poly(diffvar)

    if not chaospy.poly.caller.is_decomposed(diffvar):
        sum(differential(poly, chaospy.poly.caller.decompose(diffvar)))

    if diffvar.shape:
        return Poly([differential(poly, pol) for pol in diffvar])

    if diffvar.dim > poly.dim:
        poly = chaospy.poly.dimension.setdim(poly, diffvar.dim)
    else:
        diffvar = chaospy.poly.dimension.setdim(diffvar, poly.dim)

    qkey = diffvar.keys[0]

    core = {}
    for key in poly.keys:

        newkey = np.array(key) - np.array(qkey)

Gradient of a polynomial.

    Args:
        poly (Poly) : polynomial to take gradient of.

    Returns:
        (Poly) : The resulting gradient.

    Examples:
        >>> q0, q1, q2 = chaospy.variable(3)
        >>> poly = 2*q0 + q1*q2
        >>> print(chaospy.gradient(poly))
        [2, q2, q1]
    """
    return differential(
        poly, chaospy.poly.collection.basis(1, 1, poly.dim, sort="GR"))

>>> print(differential(poly, q1))
        [0, 0, 2q0q1]
    """
    poly = Poly(poly)
    diffvar = Poly(diffvar)

    if not chaospy.poly.caller.is_decomposed(diffvar):
        sum(differential(poly, chaospy.poly.caller.decompose(diffvar)))

    if diffvar.shape:
        return Poly([differential(poly, pol) for pol in diffvar])

    if diffvar.dim > poly.dim:
        poly = chaospy.poly.dimension.setdim(poly, diffvar.dim)
    else:
        diffvar = chaospy.poly.dimension.setdim(diffvar, poly.dim)

    qkey = diffvar.keys[0]

    core = {}
    for key in poly.keys:

        newkey = np.array(key) - np.array(qkey)

        if np.any(newkey < 0):
            continue
        newkey = tuple(newkey)
        core[newkey] = poly.A[key] * np.prod(
            [factorial(key[idx], exact=True) / factorial(newkey[idx], exact=True)
             for idx in range(poly.dim)])

    return Poly(core, poly.dim, poly.shape, poly.dtype)

Q (Poly):
            Polynomial to differentiate by. Must be decomposed. If polynomial
            array, the output is the Jacobian matrix.
    """
    P, Q = Poly(P), Poly(Q)

    if not chaospy.poly.is_decomposed(Q):
        differential(chaospy.poly.decompose(Q)).sum(0)

    if Q.shape:
        return Poly([differential(P, q) for q in Q])

    if Q.dim>P.dim:
        P = chaospy.poly.setdim(P, Q.dim)
    else:
        Q = chaospy.poly.setdim(Q, P.dim)

    qkey = Q.keys[0]

    A = {}
    for key in P.keys:

        newkey = numpy.array(key) - numpy.array(qkey)

        if numpy.any(newkey<0):
            continue

        A[tuple(newkey)] = P.A[key]*numpy.prod([factorial(key[i], \
            exact=True)/factorial(newkey[i], exact=True) \
            for i in range(P.dim)])

    return Poly(B, P.dim, P.shape, P.dtype)

def __init__(self, core=None, dim=None, shape=None, dtype=None):
        """
        Constructor for the Poly class.

        Args:
            A (numpy.ndarray, dict, Poly) : The polynomial coefficient Tensor.
                    Where A[(i,j,k)] corresponds to a_{ijk} x^i y^j z^k
                    (A[i][j][k] for list and tuple)
            dim (int) : the dimensionality of the polynomial.  Automatically
                    set if A contains a value.
            shape (tuple) : the number of polynomials represented.
                    Automatically set if A contains a value.
            dtype (type) : The type of the polynomial coefficients
        """
        core, dim, shape, dtype = chaospy.poly.constructor.preprocess(
            core, dim, shape, dtype)

        self.keys = sorted(core.keys(), key=sort_key)
        self.dim = dim
        self.shape = shape
        self.dtype = dtype
        self.A = core