@@ -1,6 +1,6 @@
 """This file is to collect those lineages that have the same condition, and have 3 or greater number of cells in their lineages."""
 
-from .import_lineage import MCF10A
+from .import_lineage import MCF10A, import_taxol
 from .LineageInputOutput import import_exp_data
 from .LineageTree import LineageTree
 from .states.StateDistributionGamma import StateDistribution as StDist
@@ -153,3 +153,15 @@ def load_replicate_groups(replicate_file_groups: list[list[str]]):
 AllLapatinib = [Lapatinib_Control + Gemcitabine_Control, Lapt25uM, Lapt50uM, Lap250uM]
 AllGemcitabine = [Lapatinib_Control + Gemcitabine_Control, Gem5uM, Gem10uM, Gem30uM]
 GFs = [pbs, egf, hgf, osm]
+
+##########################################################
+# NEW PACLITAXEL DATA
+untreated, taxol_05, taxol_1, taxol_15, taxol_2, taxol_25, taxol_3, taxol_4 = import_taxol()
+# 555, 1119, 486, 382, 272, 365, 246, 266
+untreated_taxol = [LineageTree(list_cells, E) for list_cells in untreated]  # 555 cells
+Taxol_1 = [LineageTree(list_cells, E) for list_cells in taxol_1]  # 486 cells
+Taxol_2 = [LineageTree(list_cells, E) for list_cells in taxol_2]  # 272 cells
+Taxol_3 = [LineageTree(list_cells, E) for list_cells in taxol_3]  # 246 cells
+Taxol_4 = [LineageTree(list_cells, E) for list_cells in taxol_4]  # 266 cells
+
+taxols = [untreated_taxol, Taxol_1, Taxol_2, Taxol_4]
@@ -12,33 +12,35 @@
 exe = ProcessPoolExecutor()
 
 
-def hide_observation(lineages: list, percentage: float) -> list:
+def hide_observation(lineages: list, percentage: float, rng=None) -> list:
     """Taking a list of lineages and the percentage of cells want to be masked, it marks those x% negative."""
+    rng = np.random.default_rng(rng)
     new_lineages = deepcopy(lineages)
     for new_lineage in new_lineages:
         for cell in new_lineage.output_lineage:
-            if bernoulli.rvs(p=percentage, size=1):
+            if bernoulli.rvs(p=percentage, size=1, random_state=rng):
                 # negate the cell observations to mask them
                 cell.obs = [-1 * o for o in cell.obs]
         new_lineage.obs = np.array([cell.obs for cell in new_lineage.output_lineage])
 
     return new_lineages
 
 
-def crossval(train_populations: list, num_states: np.ndarray):
+def crossval(train_populations: list, num_states: np.ndarray, rng=None):
     """Perform cross validation for the experimental data which runs in parallel for all states.
     :param train_populations: the populations after applying hide_observation. This includes the list of list of lineages.
     :param hidden_indexes: is a list of list of np.arrays for each lineage,
     filled with zeros and ones. ones refer to the index of those cells that have been hidden.
     :param hidden_obs: list of list of tuples of observations that have been masked in the train_lineage.
     :param num_states: is a range of states we want to run the cross validation for.
     """
+    rng = np.random.default_rng(rng)
     # fit training data by parallel.
 
     # save the tHMMobj for each number of states that is being run
     LLs = []
     for k in num_states:
-        tHMMobj_list, _, gamma_list = Analyze_list(train_populations, k)
+        tHMMobj_list, _, gamma_list = Analyze_list(train_populations, k, rng=rng)
 
         Logls = 0
         # calculate the log likelihood of hidden observations
 
@@ -0,0 +1,47 @@
+"""In this file we plot the raw data before fitting from paclitaxel treated HCC1143 cells"""
+
+# from ..Lineage_collections import taxols as Taxol_lin_list
+from string import ascii_lowercase
+
+from ..import_lineage import import_taxol_file
+from ..LineageTree import LineageTree
+from ..plotTree import plotLineage
+from ..states.StateDistributionGaPhs import StateDistribution
+from .common import getSetup
+
+desired_num_states = 2
+E = [StateDistribution() for _ in range(desired_num_states)]
+untreated_t = import_taxol_file()
+untreated_taxol = [LineageTree(list_cells, E) for list_cells in untreated_t]
+c = 0
+for lins in untreated_t:
+    for _cell in lins:
+        c += 1
+# counts = len(untreated_taxol)
+counts = 100
+print("number of lineages", len(untreated_taxol), " and the number of total cells: ", c)
+
+
+def makeFigure():
+    """
+    Makes figure S18.
+    """
+    titles = ["Untreated", "Taxol 1 nM", "Taxol 2 nM", "Taxol 3 nM", "Taxol 4 nM"]
+    ax, f = getSetup((17, 35), (counts, 4))
+
+    for i in range(counts):
+        ax[4 * i].axis("off")
+        # ax[4 * i + 1].axis('off')
+        # ax[4 * i + 2].axis('off')
+        # ax[4 * i + 3].axis('off')
+        plotLineage(untreated_taxol[i], ax[4 * i], color=False)
+        # plotLineage(Taxol_lin_list[1][i], ax[4 * i + 1], color=False)
+        # plotLineage(Taxol_lin_list[2][i], ax[4 * i + 2], color=False)
+        # plotLineage(Taxol_lin_list[3][i], ax[4 * i + 3], color=False)
+
+    for i in range(4):
+        ax[i].axis("off")
+        ax[i].text(-0.2, 1.55, ascii_lowercase[i], transform=ax[i].transAxes, fontsize=20, fontweight="bold", va="top")
+        ax[i].text(0.0, 1.55, titles[i], transform=ax[i].transAxes, fontsize=20, va="top")
+
+    return f
@@ -0,0 +1,37 @@
+"""BIC for Taxol data."""
+
+import pickle
+
+import numpy as np
+from matplotlib.ticker import MaxNLocator
+
+from ..Lineage_collections import taxols as Taxol_lin_list
+from .common import getSetup
+from .figure9 import find_BIC
+
+desired_num_states = np.arange(1, 8)
+
+
+def makeFigure():
+    """
+    Makes figure 9.
+    """
+    ax, f = getSetup((14, 4), (1, 3))
+
+    BIC, Obj = find_BIC(Taxol_lin_list, desired_num_states, num_cells=1041)
+
+    # create a pickle file
+    pik1 = open("taxols.pkl", "wb")
+    for taxol_tHMMobj_list in Obj:
+        for laps in taxol_tHMMobj_list:
+            pickle.dump(laps, pik1)
+    pik1.close()
+
+    i = 0
+    ax[i].set_xlabel("Number of States Predicted")
+    ax[i].set_ylabel("Normalized BIC")
+    ax[i].xaxis.set_major_locator(MaxNLocator(integer=True))
+
+    ax[0].set_title("Taxol Treated Populations")
+
+    return f