How to use the chaospy.generate_quadrature function in chaospy

def test_approx_quadrature():
    dist = cp.Iid(normal(), dim)
    nodes, weights = cp.generate_quadrature(order, dist, rule="C")

data = get_data(nrun,cdir)
#rescale oxygen flux
data[:,0,:] = -data[:,0,:]*86400.

uq = profit.UQ(yaml='uq.yaml')
distribution = cp.J(*uq.params.values())
sparse=uq.backend.sparse
if sparse:
  order=2*3
else:
  order=3+1

# actually start the postprocessing now:

nodes, weights = cp.generate_quadrature(order, distribution, rule='G',sparse=sparse)
expansion,norms = cp.orth_ttr(3, distribution,retall=True)
approx_denit = cp.fit_quadrature(expansion, nodes, weights, np.mean(data[:,1,:], axis=1))
approx_oxy = cp.fit_quadrature(expansion, nodes, weights, np.mean(data[:,0,:], axis=1))

annual_oxy = cp.fit_quadrature(expansion,nodes,weights,data[:,0,:])
annual_denit = cp.fit_quadrature(expansion,nodes,weights,data[:,1,:])

s_denit = cp.descriptives.sensitivity.Sens_m(annual_denit, distribution)
s_oxy = cp.descriptives.sensitivity.Sens_m(annual_oxy, distribution)

df_oxy = cp.Std(annual_oxy,distribution)
df_denit = cp.Std(annual_denit,distribution)
f0_oxy = cp.E(annual_oxy,distribution)
f0_denit = cp.E(annual_denit,distribution)

expected = cp.E(foo_approx, distribution)
    deviation = cp.Std(foo_approx, distribution)

    plt.fill_between(
        coord, expected-deviation, expected+deviation,
        color="k", alpha=0.3
    )
    plt.plot(coord, expected, "k--", lw=3)
    plt.xlabel(r"\verb;coord;")
    plt.ylabel(r"function evaluations \verb;foo;")
    plt.title("Results using point collocation method")
    plt.savefig("results_collocation.png")
    plt.clf()


    absissas, weights = cp.generate_quadrature(8, distribution, "C")
    evals = [foo(coord, val) for val in absissas.T]
    foo_approx = cp.fit_quadrature(polynomial_expansion, absissas, weights, evals)
    expected = cp.E(foo_approx, distribution)
    deviation = cp.Std(foo_approx, distribution)

    plt.fill_between(
        coord, expected-deviation, expected+deviation,
        color="k", alpha=0.3
    )
    plt.plot(coord, expected, "k--", lw=3)
    plt.xlabel(r"\verb;coord;")
    plt.ylabel(r"function evaluations \verb;foo;")
    plt.title("Results using psuedo-spectral projection method")
    plt.savefig("results_spectral.png")
    plt.clf()

distribution = self.create_distribution(uncertain_parameters=uncertain_parameters)


        # Create the Multivariate normal distribution
        dist_R = []
        for parameter in uncertain_parameters:
            dist_R.append(cp.Normal())

        dist_R = cp.J(*dist_R)

        P = cp.orth_ttr(polynomial_order, dist_R)

        if quadrature_order is None:
            quadrature_order = polynomial_order + 2

        nodes_R, weights_R = cp.generate_quadrature(quadrature_order,
                                                    dist_R,
                                                    rule="J",
                                                    sparse=True)


        nodes = distribution.inv(dist_R.fwd(nodes_R))
        # weights = weights_R*distribution.pdf(nodes)/dist_R.pdf(nodes_R)

        # Running the model
        data = self.runmodel.run(nodes, uncertain_parameters)

        data.method = "polynomial chaos expansion with the pseudo-spectral method and the Rosenblatt transformation. polynomial_order={}, quadrature_order={}".format(polynomial_order, quadrature_order)

        logger = get_logger(self)

        U_hat = {}

#will hold the products of PCE basis functions phi_k and lagrange interpolation polynomials a_1d
                    cross_prod = []

                    for n in range(self.N):
                        #GQ using n points can exactly integrate polynomials of order 2n-1:
                        #solve for required number of points when given order
                        n_quad_points = int(np.ceil((orders[n] + 1) / 2))

                        #generate Gaussian quad rule
                        if isinstance(self.sampler.params_distribution[n], cp.DiscreteUniform):
                            xi = self.xi_1d[n][l[n]]
                            wi = self.wi_1d[n][l[n]]
                            #TODO: remove when chaospy discrete weights have been fixed
                            # wi = np.ones(xi.size)/xi.size
                        else:
                            xi, wi = cp.generate_quadrature(n_quad_points - 1, self.sampler.params_distribution[n], rule="G")
                            xi = xi[0]

                        #number of colloc points = number of Lagrange polynomials
                        n_lagrange_poly = int(self.xi_1d[n][l[n]].size)

                        #compute the v coefficients = coefficients of SC2PCE mapping
                        v_coefs_n = []
                        for j in range(n_lagrange_poly):
                            #compute values of Lagrange polys at quadrature points
                            l_j = np.array([lagrange_poly(xi[i], self.xi_1d[n][l[n]], j) for i in range(xi.size)])
                            #each coef is the integral of the lagrange poly times the current orthogonal PCE poly
                            v_coefs_n.append(np.sum(l_j*phi_k[n](xi)*wi))
                        cross_prod.append(v_coefs_n)

                    #tensor product of all integrals
                    integrals = np.array(list(product(*cross_prod)))

else:
            wi_1d = {}

            params = self.sampler.params_distribution

            for n in range(self.N):
                # 1d weights for n-th parameter
                wi_1d[n] = {}
                # loop over all level of collocation method
                for level in range(1, self.L + 1):
                    # current SC nodes over dimension n and level
                    xi_1d = self.xi_1d[n][level]
                    wi_1d[n][level] = np.zeros(xi_1d.size)

                    # generate a quadrature rule to compute the SC weights
                    xi_quad, wi_quad = cp.generate_quadrature(level, params[n], rule=rule)
                    xi_quad = xi_quad[0]

                    # compute integral of the lagrange polynomial through xi_1d, weighted
                    # by the input distributions:
                    # w_j = int L_j(xi) p(xi) dxi j = 1,..,xi_1d.size
                    for j in range(xi_1d.size):
                        # values of L_i(xi_quad)
                        lagrange_quad = np.zeros(xi_quad.size)
                        for i in range(xi_quad.size):
                            lagrange_quad[i] = lagrange_poly(xi_quad[i], xi_1d, j)
                        # quadrature
                        wi_1d[n][level][j] = np.sum(lagrange_quad * wi_quad)

            return wi_1d

# Regression variante (Point collocation method)
        if regression:
            # Change the default rule
            if rule == "G":
                self.rule = "M"

            # Generates samples
            self.n_samples = 2 * len(self.P)
            nodes = cp.generate_samples(order=self.n_samples,
                                        domain=self.distribution,
                                        rule=self.rule)

        # Projection variante (Pseudo-spectral method)
        else:
            # Nodes and weights for the integration
            nodes, _ = cp.generate_quadrature(order=polynomial_order,
                                              dist=self.distribution,
                                              rule=self.rule,
                                              sparse=sparse,
                                              growth=self.quad_growth)
            # Number of samples
            self.n_samples = len(nodes[0])

        # Reorganize nodes according to params type: scalar (float, integer) or list
        self._nodes = []
        ipar = 0
        for j in self.params_size:
            # Scalar
            if j == 1:
                self._nodes.append(nodes[ipar:ipar + 1].flatten())
            # List
            else:

j = 0
            else:
                j = 1

            for n in range(N):
                for i in range(L):
                    xi_i, wi_i = cp.generate_quadrature(i + j,
                                                        self.params_distribution[n],
                                                        rule=self.quad_rule,
                                                        growth=self.growth)

                    self.xi_1d[n][i + 1] = xi_i[0]
                    self.wi_1d[n][i + 1] = wi_i
        else:
            for n in range(N):
                xi_i, wi_i = cp.generate_quadrature(self.polynomial_order[n],
                                                    self.params_distribution[n],
                                                    rule=self.quad_rule,
                                                    growth=self.growth)

                self.xi_1d[n][self.polynomial_order[n]] = xi_i[0]
                self.wi_1d[n][self.polynomial_order[n]] = wi_i

def evaluate_postprocessing(distribution,data,expansion):
    import matplotlib.pyplot as plt
    from profit import read_input
    from chaospy import generate_quadrature, orth_ttr, fit_quadrature, E, Std, descriptives
    
    nodes, weights = generate_quadrature(uq.backend.order+1, distribution, rule='G')
    expansion = orth_ttr(uq.backend.order, distribution)
    approx = fit_quadrature(expansion, nodes, weights, np.mean(data[:,0,:], axis=1))
    urange = list(uq.params.values())[0].range()
    vrange = list(uq.params.values())[1].range()
    u = np.linspace(urange[0], urange[1], 100)
    v = np.linspace(vrange[0], vrange[1], 100)
    U, V = np.meshgrid(u, v)
    c = approx(U,V)    

    # for 3 parameters:
    #wrange = list(uq.params.values())[2].range()
    #w = np.linspace(wrange[0], wrange[1], 100)
    #W = 0.03*np.ones(U.shape)
    #c = approx(U,V,W)
            
    plt.figure()

kws:
            Extra args passed to `chaospy.generate_quadrature`.
    """
    if order is None:
        order = int(1000./numpy.log2(len(dist)+1))
    assert isinstance(order, int)
    assert isinstance(dist, chaospy.Dist)
    k_loc = numpy.asarray(k_loc, dtype=int)
    assert k_loc.shape == (len(dist),), "incorrect size of exponents"
    assert k_loc.dtype == int, "exponents have the wrong dtype"

    if (tuple(k_loc), dist) in MOMENTS_RESULTS:
        return MOMENTS_RESULTS[tuple(k_loc), dist]

    if (tuple(k_loc), dist, order) not in MOMENTS_QUADS:
        MOMENTS_QUADS[tuple(k_loc), dist, order] = chaospy.generate_quadrature(
            order, dist, rule=rule, **kws)
    X, W = MOMENTS_QUADS[tuple(k_loc), dist, order]

    out = numpy.sum(numpy.prod(X.T**k_loc, 1)*W)
    MOMENTS_RESULTS[tuple(k_loc), dist] = out
    return out