def run_lineage(means, weights, curve_function='<STR_LIT>', curve_dimensions=<NUM_LIT:6>):
if curve_function=='<STR_LIT>':<EOL><INDENT>func = poly_curve<EOL><DEDENT>elif curve_function=='<STR_LIT>':<EOL><INDENT>func = fourier_series<EOL><DEDENT>X = dim_reduce(means, weights, <NUM_LIT:2>)<EOL>reduced_data = np.dot(X.T, weights)<EOL>if X.shape[<NUM_LIT:0>]==<NUM_LIT:2>:<EOL><INDENT>reduced_data = np.dot(X, weights)<EOL><DEDENT>cells = weights.shape[<NUM_LIT:1>]<EOL>clusters = means.shape[<NUM_LIT:1>]<EOL>cell_cluster_assignments = weights.argmax(<NUM_LIT:0>)<EOL>cluster_curves = []<EOL>cluster_fitted_vals = reduced_data.copy()<EOL>cluster_edges = []<EOL>for c in range(clusters):<EOL><INDENT>cluster_cells = reduced_data[:, cell_cluster_assignments==c]<EOL>if len(cluster_cells) == <NUM_LIT:0>:<EOL><INDENT>cluster_edges.append([])<EOL>continue<EOL><DEDENT>if cluster_cells.shape[<NUM_LIT:1>] < <NUM_LIT:2>:<EOL><INDENT>cluster_edges.append([])<EOL>continue<EOL><DEDENT>elif cluster_cells.shape[<NUM_LIT:1>] < curve_dimensions:<EOL><INDENT>tc = cluster_cells.shape[<NUM_LIT:1>]-<NUM_LIT:1><EOL><DEDENT>else:<EOL><INDENT>tc = curve_dimensions<EOL><DEDENT>if curve_function=='<STR_LIT>':<EOL><INDENT>p0 = [<NUM_LIT:1.0>]*tc<EOL>p0[<NUM_LIT:1>] = <NUM_LIT><EOL>bounds = (-np.inf, np.inf)<EOL><DEDENT>else:<EOL><INDENT>p0 = [<NUM_LIT:1.0>]*tc<EOL>bounds = (-np.inf, np.inf)<EOL><DEDENT>p_x, pcov_x = curve_fit(func, cluster_cells[<NUM_LIT:0>,:],<EOL>cluster_cells[<NUM_LIT:1>,:],<EOL>p0=p0, bounds=bounds)<EOL>perr_x = np.sum(np.sqrt(np.diag(pcov_x)))<EOL>p_y, pcov_y = curve_fit(func, cluster_cells[<NUM_LIT:1>,:],<EOL>cluster_cells[<NUM_LIT:0>,:],<EOL>p0=p0, bounds=bounds)<EOL>perr_y = np.sum(np.sqrt(np.diag(pcov_y)))<EOL>if perr_x <= perr_y:<EOL><INDENT>x_vals = reduced_data[<NUM_LIT:0>,:]<EOL>cluster_curves.append(p_x)<EOL>y_vals = np.array([func(x, *p_x) for x in x_vals])<EOL>fitted_vals = np.array([x_vals, y_vals])<EOL>cluster_fitted_vals[:,cell_cluster_assignments==c] = fitted_vals[:,cell_cluster_assignments==c]<EOL>x_indices = np.argsort(x_vals)<EOL>x_indices = [x for x in x_indices if cell_cluster_assignments[x]==c]<EOL>new_cluster_edges = []<EOL>for i, j in zip(x_indices[:-<NUM_LIT:1>], x_indices[<NUM_LIT:1>:]):<EOL><INDENT>new_cluster_edges.append((i,j))<EOL><DEDENT>cluster_edges.append(new_cluster_edges)<EOL><DEDENT>else:<EOL><INDENT>y_vals = reduced_data[<NUM_LIT:1>,:]<EOL>cluster_curves.append(p_y)<EOL>x_vals = np.array([func(x, *p_y) for x in y_vals])<EOL>fitted_vals = np.array([x_vals, y_vals])<EOL>cluster_fitted_vals[:,cell_cluster_assignments==c] = fitted_vals[:,cell_cluster_assignments==c]<EOL>y_indices = np.argsort(y_vals)<EOL>y_indices = [x for x in y_indices if cell_cluster_assignments[x]==c]<EOL>new_cluster_edges = []<EOL>for i,j in zip(y_indices[:-<NUM_LIT:1>], y_indices[<NUM_LIT:1>:]):<EOL><INDENT>new_cluster_edges.append((i,j))<EOL><DEDENT>cluster_edges.append(new_cluster_edges)<EOL><DEDENT><DEDENT>distances = squareform(pdist(cluster_fitted_vals.T))<EOL>for c1 in range(clusters):<EOL><INDENT>min_dist = np.inf<EOL>min_index = None<EOL>if sum(cell_cluster_assignments==c1)==<NUM_LIT:0>:<EOL><INDENT>continue<EOL><DEDENT>for c2 in range(clusters):<EOL><INDENT>if sum(cell_cluster_assignments==c2)==<NUM_LIT:0>:<EOL><INDENT>continue<EOL><DEDENT>if c1!=c2:<EOL><INDENT>distances_c = distances[cell_cluster_assignments==c1,:][:, cell_cluster_assignments==c2]<EOL>mindex = np.unravel_index(distances_c.argmin(), distances_c.shape)<EOL>if distances_c[mindex] < min_dist:<EOL><INDENT>min_dist = distances_c[mindex]<EOL>min_index = np.where(distances==min_dist)<EOL>min_index = (min_index[<NUM_LIT:0>][<NUM_LIT:0>], min_index[<NUM_LIT:1>][<NUM_LIT:0>])<EOL><DEDENT><DEDENT><DEDENT>cluster_edges[c1].append(min_index)<EOL><DEDENT>cluster_edges = [i for sublist in cluster_edges for i in sublist]<EOL>return cluster_curves, cluster_fitted_vals, cluster_edges, cell_cluster_assignments<EOL>
Lineage graph produced by minimum spanning tree

Args:
    means (array): genes x clusters - output of state estimation
    weights (array): clusters x cells - output of state estimation
    curve_function (string): either 'poly' or 'fourier'. Default: 'poly'
    curve_dimensions (int): number of parameters for the curve. Default: 6

Returns:
    curve parameters: list of lists for each cluster
    smoothed data in 2d space: 2 x cells
    list of edges: pairs of cell indices
    cell cluster assignments: list of ints