diff --git a/growthAnalysis/aniso.py b/growthAnalysis/aniso.py
new file mode 100644
index 0000000000000000000000000000000000000000..7aa5a620177e9109d376df3b59d63d6120f15d09
--- /dev/null
+++ b/growthAnalysis/aniso.py
@@ -0,0 +1,198 @@
+import numpy as np
+from scipy.spatial import ConvexHull, KDTree
+from scipy.linalg import lstsq
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import common.lin as lin
+import common.seg as seg
+
+def centre_points(ps):
+ return ps - np.mean(ps, axis=0)
+
+def merge_points(cell_regions, cids):
+ return np.vstack(tuple([cell_regions[cid] for cid in cids]))
+
+def fract_anisotropy(s):
+ tr = np.mean(s)
+ return (np.sqrt(1.5) * (np.sqrt(np.sum((s - tr)**2)) / np.sqrt(np.sum(s**2))))
+
+def fract_anisotropy1(s):
+ tr = np.mean(s)
+ return (np.sqrt(1.5) * (np.sqrt(np.sum((s - tr)**2))))
+
+def eival_anisotropy(s):
+ return s[0] / s[1]
+
+def hulls(regs, linss, t1, t2, i):
+ ds = linss[(t1, t2)][i]
+
+ if not ds:
+ return 0.0, np.array([]), np.array([])
+
+ ps1 = centre_points(regs[t1][i])
+ ps2 = centre_points(merge_points(regs[t2], ds))
+
+ h1 = ConvexHull(ps1)
+ h2 = ConvexHull(ps2)
+
+ vs1 = h1.points[h1.vertices]
+ vs2 = h2.points[h2.vertices]
+
+ return vs1, vs2
+
+def ganiso_f(regs, linss, t1, t2, i, descr=fract_anisotropy):
+ dt = t2 - t1
+
+ linsL = lin.filterLinssBetween(linss, (t1, t2))
+
+ ds = lin.tSearchN(linsL, i, len(linsL), 0)
+
+ if not ds:
+ return None, np.array([]), np.array([])
+
+ ps1 = centre_points(regs[t1][i])
+ ps2 = centre_points(merge_points(regs[t2], ds))
+
+ h1 = ConvexHull(ps1)
+ h2 = ConvexHull(ps2)
+
+ vs1 = h1.points[h1.vertices]
+ vs2 = h2.points[h2.vertices]
+
+ vs2_q = KDTree(vs2)
+ _, ids = vs2_q.query(vs1)
+ vs2_ = vs2[ids, :]
+
+ A = lstsq(vs1, vs2_)
+
+ _, s, _ = np.linalg.svd(A[0])
+
+ return (descr(s) / dt), vs1, vs2
+
+def get_eivals(D):
+ cv = np.cov(D)
+ eivals, eivecs = np.linalg.eigh(cv)
+ key = np.argsort(eivals)[::-1][:]
+ eivals, eivecs = eivals[key], eivecs[:, key]
+
+ return eivals, eivecs
+
+def sample_points(D):
+ return D[0:len(D):2, :]
+
+def ganiso_f_pca(regs, linss, t1, t2, i):
+ dt = t2 - t1
+
+ ds = linss[(t1, t2)][i]
+
+ if not ds:
+ return 0.0, np.array([]), np.array([])
+
+ ps1 = sample_points(centre_points(regs[t1][i]))
+ ps2 = sample_points(centre_points(merge_points(regs[t2], ds)))
+
+ eivals1, _ = get_eivals(ps1.T)
+ eivals2, _ = get_eivals(ps2.T)
+
+ a1 = eival_anisotropy(eivals1)
+ a2 = eival_anisotropy(eivals2)
+
+ return (((a2-a1) / a1) / dt), ps1, ps2
+
+def plot_points3(vs1, vs2):
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+
+ ax.scatter(vs1[:, 0], vs1[:, 1], vs1[:, 2], c='b')
+ ax.scatter(vs2[:, 0], vs2[:, 1], vs2[:, 2], c='r')
+
+ plt.show()
+
+def plot_points2(vs1, vs2):
+ fig = plt.figure()
+ ax = fig.add_subplot(111)
+
+ ax.scatter(vs1[:, 0], vs1[:, 1], c='b')
+ ax.scatter(vs3[:, 0], vs3[:, 1], c='g')
+
+ for i in range(len(vs1)):
+ ax.annotate(str(i), (vs1[i, 0], vs1[i, 1]))
+
+ for i in range(len(vs2)):
+ ax.annotate(str(i), (vs2[i, 0], vs2[i, 1]))
+
+ plt.show()
+
+def ganisos_forward(regs, tss, linss, t1, t2, descr=fract_anisotropy):
+ ans_ = []
+ for c in tss[t1]:
+ an, _, _ = ganiso_f(regs, linss, t1, t2, c.cid, descr=descr)
+ if an:
+ ans_.append((c.cid, an))
+
+ return dict(ans_)
+
+def ganisos_backward(regs, tss, linss, t1, t2):
+ linsL = lin.filterLinssBetween(linss, (t1, t2))
+ ans_ = []
+
+ for c in tss[t1]:
+ an, _, _ = ganiso_f(regs, linss, t1, t2, c.cid)
+ ds = lin.tSearchN(linsL, c.cid, len(linsL), 0)
+ for d in ds:
+ ans_.append((d, an))
+
+ return dict(ans_)
+
+def ganisos_backward_data(gans, tss, linss, t1, t2):
+ ans_ = []
+ im = lin.invertMap(linss[(t1, t2)])
+ for d in tss[t2].cells.keys():
+ i = im.get(d, -1)
+ if i > 0:
+ an = gans[t1].get(i, 0.0)
+ ans_.append((d, an))
+ else:
+ ans_.append((d, 0.0))
+
+ return dict(ans_)
+
+def plot_cells_aniso(ts, ans):
+ ps = np.array([[c.pos.x, c.pos.y, c.pos.z] for c in ts])
+ cs = [ans.get(c.cid, 0.0) for c in ts]
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+
+ ax.scatter(ps[:, 0], ps[:, 1], ps[:, 2], c=cs, s=100, linewidths=0.5,
+ cmap=plt.cm.hot, alpha=0.9, edgecolors='black', vmin=0.1, vmax=0.8)
+
+ plt.show()
+
+
+def plot_pca(D):
+ Dc = centre_points(D)
+ eivals, eivecs = get_eivals(Dc.T)
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+
+ ax.scatter(Dc[:, 0], Dc[:, 1], Dc[:, 2], c='b')
+ ax.plot([0, eivals[0]*eivecs[0, 0]],
+ [0, eivals[0]*eivecs[0, 1]],
+ [0, eivals[0]*eivecs[0, 2]],
+ c='r')
+
+ ax.plot([0, eivals[1]*eivecs[1, 0]],
+ [0, eivals[1]*eivecs[1, 1]],
+ [0, eivals[1]*eivecs[1, 2]],
+ c='r')
+
+ ax.plot([0, eivals[2]*eivecs[2, 0]],
+ [0, eivals[2]*eivecs[2, 1]],
+ [0, eivals[2]*eivecs[2, 2]],
+ c='r')
+
+ plt.show()
+
+ return eivecs, eivals
diff --git a/growthAnalysis/anisos.ipynb b/growthAnalysis/anisos.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..e55b858f380930deb202df8bb9bd5c3aeb945f97
--- /dev/null
+++ b/growthAnalysis/anisos.ipynb
@@ -0,0 +1,405 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import patternTransitions as p\n",
+ "from itertools import repeat\n",
+ "import seaborn as sns\n",
+ "import matplotlib.pyplot as plt\n",
+ "from mpl_toolkits.mplot3d import Axes3D\n",
+ "import common.seg as seg\n",
+ "import _pickle as cPickle\n",
+ "import aniso\n",
+ "import numpy as np\n",
+ "import common.lin as lin\n",
+ "import matplotlib.ticker as ticker"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_params = {'01' : {'elev': 120, 'azim' : 270, 'vmax':0.4},\n",
+ " '02' : {'elev': 70, 'azim' : -20, 'vmax':0.3},\n",
+ " '03' : {'elev': 90, 'azim' : 270, 'vmax':0.2},\n",
+ " '04' : {'elev':80, 'azim': 50, 'vmax': 0.2},\n",
+ " '05' : {'elev':70, 'azim': 220, 'vmax':0.2},\n",
+ " '06' : {'elev':90, 'azim': 270, 'vmax':0.2}}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def mkViolinPlot(d, ax, ws=0.7):\n",
+ " xs_labs = sorted(d.keys())\n",
+ " xs = range(len(xs_labs))\n",
+ " ys = [d[x] for x in xs_labs]\n",
+ " \n",
+ " ax.yaxis.set_major_locator(ticker.MultipleLocator(0.005))\n",
+ " sns.set_style(\"whitegrid\", {'grid.linestyle':'--'})\n",
+ " cls = sns.color_palette()\n",
+ " ax.set_xticks(xs)\n",
+ " ax.set_xticklabels(xs_labs)\n",
+ " violin_parts = ax.violinplot(ys, xs, showmeans=True, widths=ws)\n",
+ " for pc in violin_parts['bodies']:\n",
+ " pc.set_facecolor(\"#f5d742\")\n",
+ " pc.set_edgecolor('black')\n",
+ "\n",
+ " for partname in ('cbars','cmins','cmaxes','cmeans'):\n",
+ " vp = violin_parts[partname]\n",
+ " vp.set_edgecolor('red')\n",
+ " vp.set_linewidth(0.5)\n",
+ " \n",
+ " plt.xticks(rotation=90)\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#fid='01'\n",
+ "#tssFn = \"lins_allflowers/f\" + fid + \"_tss.pkl\"\n",
+ "#linssFn = \"lins_allflowers/f\" + fid + \"_linss.pkl\"\n",
+ "\n",
+ "#with open(tssFn, 'rb') as tssF:\n",
+ "# tss = cPickle.load(tssF, encoding='latin1')\n",
+ "\n",
+ "#with open(linssFn, 'rb') as linssF:\n",
+ "# linss = cPickle.load(linssF, encoding='latin1')\n",
+ "\n",
+ "#if fid=='01':\n",
+ "# lin.filterL1(tss)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dataDir = \"../data/\"\n",
+ "tss, linss = lin.mkSeries1(dataDir)\n",
+ "lin.filterL1(tss)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import importlib\n",
+ "importlib.reload(seg)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Per pattern"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fid='01'\n",
+ "G = p.mkTGraphN()\n",
+ "regs = seg.getCellRegions(dataDir=dataDir, fid=fid, d=3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gans = p.calcAnisotropies(regs, tss, linss)\n",
+ "gans_pat = p.addCombPatterns(gans, tss)\n",
+ "gans_state = p.getGAnisosPerPattern(G, gans_pat)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.rcParams.update({'font.size': 12})\n",
+ "fig = plt.figure(figsize=(9, 4.5))\n",
+ "ax = fig.add_subplot('111')\n",
+ "\n",
+ "ax.set_xlabel(\"cell state\")\n",
+ "ax.set_ylabel(r'anisotropy rate ($h^{-1}$)')\n",
+ "ax.yaxis.set_major_locator(ticker.MultipleLocator(0.05))\n",
+ "\n",
+ "mkViolinPlot(gans_state, ax)\n",
+ "\n",
+ "fout = \"pattern_anisotropy.png\"\n",
+ "plt.savefig(fout, dpi=300, bbox_inches='tight')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def tToStage(t):\n",
+ " if t==10: return 0\n",
+ " elif t==40: return 1\n",
+ " elif t==96: return 2\n",
+ " elif t==120: return 3\n",
+ " elif t==132: return 4"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Per pattern on graph"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "G = p.mkTGraphN()\n",
+ "p.addGAnisos(G, gans_pat)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p.drawTGraph(G, p.ganiso)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Per gene"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gans_gene = p.getGAnisosPerGene(tss, gans)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig = plt.figure(figsize=(10, 5))\n",
+ "ax = fig.add_subplot('111')\n",
+ "\n",
+ "ax.set_xlabel(\"gene\")\n",
+ "ax.set_ylabel(r'anisotropy ($h^{-1}$)')\n",
+ "ax.yaxis.set_major_locator(ticker.MultipleLocator(0.05))\n",
+ "\n",
+ "mkViolinPlot(gans_gene, ax)\n",
+ "\n",
+ "fout = \"gene_anisotropy.png\"\n",
+ "plt.savefig(fout, dpi=300, bbox_inches='tight')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Per gene, per stage"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gans_gene_time = p.getGAnisosPerGeneTime(tss, gans)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for t in sorted(gans_gene_time.keys()):\n",
+ " fig = plt.figure(figsize=(10, 5))\n",
+ " ax = fig.add_subplot('111')\n",
+ " \n",
+ " ax.set_title(\"stage {t}\".format(t=tToStage(t)))\n",
+ " ax.set_xlabel(\"gene\")\n",
+ " ax.set_ylabel(r'anisotropy ($h^{-1}$)')\n",
+ " ax.set_ylim((0.0, 0.03))\n",
+ " ax.yaxis.set_major_locator(ticker.MultipleLocator(0.05))\n",
+ "\n",
+ " mkViolinPlot(gans_gene_time[t], ax, ws=0.4)\n",
+ "\n",
+ " fout = \"gene_anisotropy_{t}.png\".format(t=t)\n",
+ " plt.savefig(fout, dpi=300, bbox_inches='tight')\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Anisotropy distributions (all to all timepoints)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fid='01'\n",
+ "#regs = seg.getCellRegions(fid=fid, d=3)\n",
+ "#tss_all = cPickle.load(file(\"lins_allflowers/f\" + fid + \"_tss.pkl\"))\n",
+ "#linss = cPickle.load(file(\"lins_allflowers/f\" + fid + \"_linss.pkl\"))\n",
+ "\n",
+ "if fid == '01':\n",
+ " lin.filterL1(tss_all)\n",
+ " tss = tss_all\n",
+ "else:\n",
+ " tss = dict([(t, seg.filterL1(ts)) for t, ts in tss_all.iteritems()])\n",
+ " \n",
+ "timepoints = sorted(tss.keys())[::2]\n",
+ "n = len(timepoints)\n",
+ "\n",
+ "fig = plt.figure(figsize=(10, 8))\n",
+ "plt.subplots_adjust(hspace=0.85, wspace=0.5)\n",
+ "\n",
+ "for i, t1 in enumerate(timepoints):\n",
+ " for k, t2 in enumerate(timepoints):\n",
+ " print(t1, end=\",\")\n",
+ " print(t2, end=\" \")\n",
+ " if t1 > t2:\n",
+ " ax = plt.subplot2grid((n, n), (i, k))\n",
+ " #plot forward\n",
+ " plt.rcParams.update({'font.size': 8})\n",
+ " ax.set_title(str(t2) + \"h\" + \"-\" + str(t1)+\"h\")\n",
+ " plt.rcParams.update({'font.size': 10})\n",
+ " ax.set_xlim([-0.05, 1.05])\n",
+ " ax.set_xticks([0.0, 0.5, 1.0])\n",
+ " ans = aniso.ganisos_forward(regs, tss, linss, t2, t1)\n",
+ " ans_ = np.array([an for _, an in ans.items() if an])\n",
+ " sns.distplot(ans_, ax=ax) \n",
+ " \n",
+ "plt.savefig(\"aniso_calcs/anisotropy_rates_f\" + str(fid) + \"_allTs.png\", dpi=300, bbox_inches='tight')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Anisotropy distributions (all to all timepoints) on templates"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_ans(ts, ans, ax, s=30, vmax=0.5):\n",
+ " xs = [c.pos.x for c in ts]\n",
+ " ys = [c.pos.y for c in ts]\n",
+ " zs = [c.pos.z for c in ts]\n",
+ " cs = [ans.get(c.cid, 0.0) for c in ts]\n",
+ " vmax = max(cs)\n",
+ " \n",
+ " \n",
+ " points = ax.scatter(xs, ys, zs, c=cs, cmap=plt.cm.jet, alpha=0.9, s=s, vmin=0.0, vmax=vmax,\n",
+ " edgecolors='black', linewidths=0.5)\n",
+ " plt.colorbar(points, shrink=0.5, ticks=[0.0, vmax])\n",
+ " plt.axis('off')\n",
+ " \n",
+ " return\n",
+ "\n",
+ "fig = plt.figure(figsize=(20, 15))\n",
+ "timepoints = sorted(tss.keys())\n",
+ "n=len(timepoints)\n",
+ "\n",
+ "for i, t1 in enumerate(timepoints):\n",
+ " for k, t2 in enumerate(timepoints):\n",
+ " print(t1, t2, end=\" \")\n",
+ " if t1 < t2:\n",
+ " ax = plt.subplot2grid((n, n), (i, k), projection='3d')\n",
+ " ax.set_title(str(t1)+\"h\" + \"-\" + str(t2) + \"h\")\n",
+ " ax.view_init(elev=plot_params[fid]['elev'], azim=plot_params[fid]['azim'])\n",
+ " #plot\n",
+ " ans = aniso.ganisos_forward(regs, tss, linss, t1, t2)\n",
+ " ans_ = dict([(cid, an) for cid, an in ans.items() if an])\n",
+ " plot_ans(tss[t1], ans_, ax, s=60)\n",
+ " elif t1 > t2:\n",
+ " ax = plt.subplot2grid((n, n), (i, k), projection='3d')\n",
+ " ax.set_title(str(t2)+\"h\" + \"-\" + str(t1) + \"h\")\n",
+ " ax.view_init(elev=plot_params[fid]['elev'], azim=plot_params[fid]['azim'])\n",
+ " #plot\n",
+ " ans = aniso.ganisos_backward(regs, tss, linss, t2, t1)\n",
+ " ans_ = dict([(cid, an) for cid, an in ans.items() if an])\n",
+ " plot_ans(tss[t1], ans_, ax, s=60)\n",
+ " else:\n",
+ " ax = plt.subplot2grid((n, n), (i, k))\n",
+ " ax.set_xlim(0.0, 1.0)\n",
+ " ax.set_ylim(0.0, 1.0)\n",
+ " ax.text(0.4, 0.5, str(t1)+\"h\", transform = ax.transAxes, fontsize=9)\n",
+ " ax.axis('off')\n",
+ " \n",
+ "fout = \"anisotropy_allTs_{fid}.png\".format(fid=fid)\n",
+ "plt.savefig(fout, dpi=300, bbox_inches='tight')"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/growthAnalysis/grates.py b/growthAnalysis/grates.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b87d7ff4781e4f01bdb17fd4f3e0bd10963dc32
--- /dev/null
+++ b/growthAnalysis/grates.py
@@ -0,0 +1,77 @@
+import numpy as np
+import common.lin as lin
+import common.edict as edict
+
+def fst(x): return x[0]
+
+def snd(x): return x[1]
+
+def avg2(v1, v2):
+ return np.mean([v1, v2])
+
+def calcVolumes(tss, linss, t1, t2):
+ tsL = lin.filterTssBetween(tss, (t1, t2))
+ linsL = lin.filterLinssBetween(linss, (t1, t2))
+ volss = lin.mergeCellVolsRefN(tsL, linsL)
+
+ cids = tsL[0].cells.keys()
+ cids = [cid for cid in cids if not volss[-1][cid] == 0]
+
+ volss1 = dict([(cid, [vols[cid] for vols in volss]) for cid in cids])
+
+ return volss1
+
+def grates_forward(tss, linss, t1, t2):
+ volss1 = calcVolumes(tss, linss, t1, t2)
+
+ grs = dict()
+ for c, vols_c in volss1.items():
+ grs[c] = ((np.log(vols_c[-1])-np.log(vols_c[0]))) / (t2-t1)
+
+ return grs
+
+def grates_backward(tss, linss, t1, t2):
+ linsL = lin.filterLinssBetween(linss, (t1, t2))
+ volss1 = calcVolumes(tss, linss, t1, t2)
+
+ grs = dict()
+ for c, vols_c in volss1.items():
+ ds = lin.tSearchN(linsL, c, len(linsL), 0)
+ for d in ds:
+ grs[d] = ((np.log(vols_c[-1])-np.log(vols_c[0]))) / (t2-t1)
+
+ return grs
+
+def grates_forward_cons(tss, linss):
+ tpoints = sorted(tss.keys())
+ grs_f = dict([(t, grates_forward(tss, linss, t, lin.succ(t)))
+ for t in tpoints if lin.succ(t)])
+
+ return grs_f
+
+def grates_backward_cons(tss, linss):
+ tpoints = sorted(tss.keys())
+ grs_b = dict([(t, grates_backward(tss, linss, lin.prev(t), t))
+ for t in tpoints if lin.prev(t)])
+
+ return grs_b
+
+def grates_avg_cons(tss, linss):
+ tpoints = sorted(tss.keys())
+ grs_f = grates_forward_cons(tss, linss)
+ grs_b = grates_backward_cons(tss, linss)
+
+ grs_avg = edict.unionWith(grs_f, grs_b,
+ lambda m1, m2: edict.unionWith(m1, m2, avg2))
+
+ return grs_avg
+
+def write_grs_mnet(t, d):
+ header = "type:float"
+ content = "\n".join([header] +
+ ["{t},{k}:{v}".format(t=t, k=k, v=v)
+ for k, v in d.items()])
+
+ return content
+
+
diff --git a/growthAnalysis/patternTransitions.py b/growthAnalysis/patternTransitions.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a8ed09e8604398cfd38f9a328b79cdc6ed821f2
--- /dev/null
+++ b/growthAnalysis/patternTransitions.py
@@ -0,0 +1,549 @@
+import _pickle as cPickle
+import networkx as nx
+import matplotlib.pyplot as plt
+import numpy as np
+import aniso as ans
+import common.seg as seg
+import common.lin as lin
+import statesN as f
+from collections import defaultdict
+import common.edict as edict
+
+patsToClusters = {2 : 1,
+ 3 : 1,
+ 6 : 1,
+ 7 : 1,
+ 10 : 1,
+ 18 : 2,
+ 27 : 2,
+ 1 : 3,
+ 32 : 3,
+ 5 : 3,
+ 9 : 4,
+ 17: 5,
+ 29 : 5,
+ 16 : 5,
+ 28 :5,
+ 11 : 6,
+ 8 : 6,
+ 12 : 7,
+ 21 :8,
+ 31:8,
+ 20:8,
+ 30:9,
+ 19:10,
+ 25:11,
+ 13:11,
+ 14:12,
+ 15:12,
+ 22:13,
+ 26:14,
+ 23:14,
+ 24:14}
+
+
+
+def findPattern1(cid, t, tss):
+ c = tss[t].cells.get(cid, None)
+
+ if c:
+ return f.get_state(c.getOnGenes())
+ else:
+ return -1
+
+def getGAnisos():
+ import gzip
+ fname = "../data/data_Jonathan/YR_01_ATLAS_iso-cell_features2.pkz"
+ f = gzip.open(fname, 'r')
+ feat_dict = cPickle.load(f)
+ f.close()
+
+ growth_anisotropy_3D_120_to_132 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 3:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_120_to_132[cid] = v
+
+ growth_anisotropy_3D_96_to_120 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 2:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_96_to_120[cid] = v
+
+
+ growth_anisotropy_3D_40_to_96 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 1:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_40_to_96[cid] = v
+
+ growth_anisotropy_3D_10_to_40 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 0:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_10_to_40[cid] = v
+
+
+ return {10 : addPats1(growth_anisotropy_3D_10_to_40, 10),
+ 40 : addPats1(growth_anisotropy_3D_40_to_96, 40),
+ 96 : addPats1(growth_anisotropy_3D_96_to_120, 96),
+ 120: addPats1(growth_anisotropy_3D_120_to_132, 120)}
+
+def getGAnisosData():
+ import gzip
+ fname = "../data/data_Jonathan/YR_01_ATLAS_iso-cell_features2.pkz"
+ f = gzip.open(fname, 'r')
+ feat_dict = cPickle.load(f)
+ f.close()
+
+ growth_anisotropy_3D_120_to_132 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 3:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_120_to_132[cid] = v
+
+ growth_anisotropy_3D_96_to_120 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 2:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_96_to_120[cid] = v
+
+
+ growth_anisotropy_3D_40_to_96 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 1:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_40_to_96[cid] = v
+
+ growth_anisotropy_3D_10_to_40 = dict()
+ for k, v in feat_dict["growth_anisotropy_3D"].iteritems():
+ if feat_dict["tp_index"][k] == 0:
+ cid = feat_dict["vid2label"][k]
+ growth_anisotropy_3D_10_to_40[cid] = v
+
+ return {10 : growth_anisotropy_3D_10_to_40,
+ 40 : growth_anisotropy_3D_40_to_96,
+ 96 : growth_anisotropy_3D_96_to_120,
+ 120: growth_anisotropy_3D_120_to_132}
+
+def calcAnisotropies(regs, tss, linss):
+ tpoints = sorted(tss.keys())
+
+ gansF = dict([(t, ans.ganisos_forward(regs, tss, linss, t, lin.succ(t)))
+ for t in tpoints if lin.succ(t)])
+
+ gansB = dict([(t, ans.ganisos_backward(regs, tss, linss, lin.prev(t), t))
+ for t in tpoints if lin.prev(t)])
+
+ gansAvg = edict.unionWith(gansF, gansB,
+ lambda m1, m2: edict.unionWith(m1, m2, avg2))
+
+ return gansAvg
+
+def addCombPatterns(gans, tss):
+ return edict.mapWithKey(gans, lambda t, ans: addPats1(ans, t, tss))
+
+def getGAnisos1():
+ tpoints = [10, 40, 96, 120]
+ regs = seg.getCellRegions(d=3)
+ tss, linss = lin.mkSeries()
+
+ gans = dict([(t, ans.ganisos_forward(regs, tss, linss, t, lin.succ(t)))
+ for t in tpoints])
+
+ return dict([(t, addPats1(gans[t], t)) for t in tpoints])
+
+def getGAnisosRev1():
+ tpoints = [40, 96, 120, 132]
+ regs = seg.getCellRegions(d=3)
+ tss, linss = lin.mkSeries()
+
+ gans = dict([(t, ans.ganisos_backward(regs, tss, linss, lin.prev(t), t))
+ for t in tpoints])
+
+ return dict([(t, addPats1(gans[t], t)) for t in tpoints])
+
+def getGAnisosRevData():
+ tpoints = [40, 96, 120, 132]
+ tss, linss = lin.mkSeries()
+ gansDataF = getGAnisosData()
+
+ gans = dict([(t, ans.ganisos_backward_data(gansDataF, tss, linss, lin.prev(t), t))
+ for t in tpoints])
+
+ return dict([(t, addPats1(gans[t], t)) for t in tpoints])
+
+
+class TPat(object):
+ def __init__(self, pid, t):
+ self.pid = pid
+ self.t = t
+ self.gr = 0
+ self.ga = 0
+
+ def addGrowthRate(self, gr):
+ self.gr = gr
+ return
+
+ def addGAniso(self, ga):
+ self.ga = ga
+ return
+
+ def __repr__(self):
+ return "(" + str(self.pid) + "," + str(self.t) + ")"
+
+ def __eq__(self, other):
+ return self.pid == other.pid and self.t == other.t
+
+ def __hash__(self):
+ return self.pid + self.t
+
+
+def toTLabel(t):
+ return str(t) + "h"
+
+def addPats(finRates, t):
+ fin = open(finRates, 'r')
+ crates = cPickle.load(fin)
+
+ return [(cid, rate, findPattern1(cid, toTLabel(t)))
+ for cid, rate in crates.iteritems()]
+
+def addPats1(cs, t, tss):
+ return [(cid, attr, findPattern1(cid, t, tss))
+ for cid, attr in cs.items()]
+
+def summRates(ratesPats, patId, f):
+ return f([rate for cid, rate, pid in ratesPats if pid == patId and rate])
+
+def gatherRates(ratesPats, patId):
+ return [rate for cid, rate, pid in ratesPats if pid == patId and rate]
+
+
+def mkPRates(ts):
+ rateFs = {10: "../data/growthRates/10h_to_40h_growthrates.pkl",
+ 40: "../data/growthRates/40h_to_96h_growthrates.pkl",
+ 96 : "../data/growthRates/96h_to_120h_growthrates.pkl",
+ 120 : "../data/growthRates/120h_to_132h_growthrates.pkl"}
+
+ return dict([(t, addPats(rateFs[t], t)) for t in ts])
+
+def mkPRates1(tss):
+ return dict([(t, addPats1(cs, t)) for t, cs in tss.iteritems()])
+
+def mkPRatesRev(ts):
+ rateFs = {40: "../data/growthRates/40h_to_10h_growthrates.pkl",
+ 96: "../data/growthRates/96h_to_40h_growthrates.pkl",
+ 120 : "../data/growthRates/120h_to_96h_growthrates.pkl",
+ 132 : "../data/growthRates/132h_to_120h_growthrates.pkl"}
+
+ return dict([(t, addPats(rateFs[t], t)) for t in ts])
+
+def addGrowthRates(G):
+ prates = mkPRates([10, 40, 96, 120])
+
+ for n in G.nodes():
+ if not n.t == 132:
+ n.addGrowthRate(summRates(prates[n.t], n.pid, np.mean))
+
+ return
+
+def addGrowthRates1(G, prates):
+ for n in G.nodes():
+ n.addGrowthRate(summRates(prates[n.t], n.pid, np.mean))
+
+ return
+
+def addGrowthRatesRev(G):
+ prates = mkPRatesRev([40, 96, 120, 132])
+
+ for n in G.nodes():
+ if not n.t == 10:
+ n.addGrowthRate(summRates(prates[n.t], n.pid, avgOr))
+
+ return
+
+def addGrowthRatesRevAvg(G):
+ prates = mkPRates([10, 40, 96, 120])
+ pratesRev = mkPRatesRev([40, 96, 120, 132])
+
+ for n in G.nodes():
+ if n.t == 10:
+ n.addGrowthRate(summRates(prates[n.t], n.pid, np.mean))
+ elif n.t == 132:
+ n.addGrowthRate(summRates(pratesRev[n.t], n.pid, avgOr))
+ else:
+ gr = summRates(prates[n.t], n.pid, np.mean)
+ grRev = summRates(pratesRev[n.t], n.pid, avgOr)
+ n.addGrowthRate(np.mean([gr, grRev]))
+
+ return
+
+def addGAnisos(G, ganisos):
+ for n in G.nodes():
+ n.addGAniso(summRates(ganisos.get(n.t, []), n.pid, avgOr))
+
+ return
+
+def addGRates(G, grates):
+ for n in G.nodes():
+ n.addGrowthRate(summRates(grates.get(n.t, []), n.pid, avgOr))
+
+ return
+
+def addGAnisos1(G):
+ ganisos = getGAnisos1()
+
+ for n in G.nodes():
+ if not n.t == 132:
+ n.addGAniso(summRates(ganisos[n.t], n.pid, avgOr))
+
+ return
+
+def addGAnisosRevAvg(G):
+ ganisos = getGAnisos1()
+ ganisosRev = getGAnisosRev1()
+
+ for n in G.nodes():
+ if n.t == 10:
+ n.addGAniso(summRates(ganisos[n.t], n.pid, np.mean))
+ elif n.t == 132:
+ n.addGAniso(summRates(ganisosRev[n.t], n.pid, avgOr))
+ else:
+ gan = summRates(ganisos[n.t], n.pid, np.mean)
+ ganRev = summRates(ganisosRev[n.t], n.pid, avgOr)
+ n.addGAniso(np.mean([gan, ganRev]))
+
+ return
+
+def getGAnisosPerPattern(G, ganisos):
+ ganisosState = defaultdict(list)
+
+ for n in G.nodes():
+ ganisosState[n.pid] += gatherRates(ganisos.get(n.t, []), n.pid)
+
+ return ganisosState
+
+def getGAnisosPerPatternTime(G, ganisos):
+ ganisosStateT = defaultdict(dict)
+
+ for n in G.nodes():
+ ganisosStateT[n.t][n.pid] = gatherRates(ganisos.get(n.t, []), n.pid)
+
+ return ganisosStateT
+
+def getGAnisosPerGene(tss, ganisos):
+ ganisosGene = defaultdict(list)
+
+ for t, gas in ganisos.items():
+ for c, a in gas.items():
+ if c in tss[t].cells:
+ for g, expr in tss[t].cells[c].exprs.items():
+ if expr: ganisosGene[g].append(a)
+
+ return ganisosGene
+
+def getGAnisosPerGeneTime(tss, ganisos):
+ ganisosGeneT = dict()
+ for t in ganisos.keys():
+ ganisosGeneT[t] = defaultdict(list)
+
+ for t, gas in ganisos.items():
+ for c, a in gas.items():
+ if c in tss[t].cells:
+ for g, expr in tss[t].cells[c].exprs.items():
+ if expr: ganisosGeneT[t][g].append(a)
+
+ return ganisosGeneT
+
+def addGAnisosRevAvgData(G):
+ ganisos = getGAnisos()
+ ganisosRev = getGAnisosRevData()
+
+ for n in G.nodes():
+ if n.t == 10:
+ n.addGAniso(summRates(ganisos[n.t], n.pid, np.mean))
+ elif n.t == 132:
+ n.addGAniso(summRates(ganisosRev[n.t], n.pid, avgOr))
+ else:
+ gan = summRates(ganisos[n.t], n.pid, np.mean)
+ ganRev = summRates(ganisosRev[n.t], n.pid, avgOr)
+ n.addGAniso(np.mean([gan, ganRev]))
+
+ return
+
+
+#blue:1
+#red:2
+#dot:3
+def mkTGraphN():
+ edges = [(TPat(1, 10), TPat(5, 40), 1),
+ (TPat(2, 10), TPat(7, 40), 1),
+ (TPat(3, 10), TPat(7, 40), 1),
+ (TPat(4, 10), TPat(4, 40), 1),
+ (TPat(31, 10), TPat(5, 40), 2),
+ (TPat(2, 10), TPat(5, 40), 2),
+ (TPat(4, 40), TPat(12, 96), 1),
+ (TPat(7, 40), TPat(12, 96), 3),
+ (TPat(7, 40), TPat(8, 96), 3),
+ (TPat(7, 40), TPat(9, 96), 1),
+ (TPat(7, 40), TPat(11, 96), 3),
+ (TPat(7, 40), TPat(10, 96), 1),
+ (TPat(5, 40), TPat(9, 96), 2),
+ (TPat(11, 96), TPat(17, 120), 1),
+ (TPat(9, 96), TPat(17, 120), 1),
+ (TPat(9, 96), TPat(16, 120), 1),
+ (TPat(8, 96), TPat(16, 120), 1),
+ (TPat(9, 96), TPat(20, 120), 3),
+ (TPat(8, 96), TPat(18, 120), 1),
+ (TPat(12, 96), TPat(18, 120), 3),
+ (TPat(12, 96), TPat(20, 120), 3),
+ (TPat(12, 96), TPat(19, 120), 1),
+ (TPat(10, 96), TPat(20, 120), 3),
+ (TPat(10, 96), TPat(19, 120), 1),
+ (TPat(10, 96), TPat(14, 120), 3),
+ (TPat(10, 96), TPat(15, 120), 1),
+ (TPat(10, 96), TPat(13, 120), 1),
+ (TPat(8, 96), TPat(17, 120), 2),
+ (TPat(13, 120), TPat(22, 132), 1),
+ (TPat(13, 120), TPat(24, 132), 1),
+ (TPat(15, 120), TPat(24, 132), 1),
+ (TPat(15, 120), TPat(23, 132), 1),
+ (TPat(14, 120), TPat(14, 132), 1),
+ (TPat(16, 120), TPat(16, 132), 3),
+ (TPat(16, 120), TPat(26, 132), 3),
+ (TPat(16, 120), TPat(27, 132), 1),
+ (TPat(17, 120), TPat(27, 132), 1),
+ (TPat(17, 120), TPat(18, 132), 1),
+ (TPat(18, 120), TPat(18, 132), 1),
+ (TPat(18, 120), TPat(26, 132), 3),
+ (TPat(19, 120), TPat(25, 132), 1),
+ (TPat(19, 120), TPat(29, 132), 1),
+ (TPat(20, 120), TPat(30, 132), 1)]
+
+ G = nx.DiGraph()
+ for n1, n2, connTyp in edges:
+ G.add_edge(n1, n2, connTyp=connTyp)
+
+ return G
+
+def nodePoss(G):
+ poss = [((1, 10), (0, -10)),
+ ((2, 10), (0, -20)),
+ ((3, 10), (0, -30)),
+ ((4, 10), (0, 0)),
+ ((32, 10), (0, -5)),
+ ((4, 40), (20, 5)),
+ ((5, 40), (20, -5)),
+ ((7, 40), (20, -18)),
+ ((8, 96), (40, -25)),
+ ((9, 96), (40, -5)),
+ ((10, 96), (40, -35)),
+ ((11, 96), (40, -15)),
+ ((12, 96), (40, 5)),
+ ((13, 120), (60, -35)),
+ ((14, 120), (60, -55)),
+ ((15, 120), (60, -45)),
+ ((16, 120), (60, -15)),
+ ((17, 120), (60, 5)),
+ ((18, 120), (60, -5)),
+ ((19, 120), (60, -25)),
+ ((20, 120), (60, 15)),
+ ((16, 132), (80, -15)),
+ ((29, 132), (80, -25)),
+ ((18, 132), (80, 15)),
+ ((22, 132), (80, -45)),
+ ((23, 132), (80, -65)),
+ ((24, 132), (80, -55)),
+ ((25, 132), (80, -35)),
+ ((26, 132), (80, -5)),
+ ((27, 132), (80, 5)),
+ ((30, 132), (80, 25)),
+ ((14, 132), (80, -75)),
+ ((31, 10), (0, 10))]
+
+ possD = dict([(tp, np.array(pos)) for tp, pos in poss])
+
+ return dict([(g, possD[(g.pid, g.t)]) for g in G.nodes()])
+
+
+def nodePoss1(G):
+ poss = [((1, 10), (0, -20)),
+ ((2, 10), (0, -40)),
+ ((3, 10), (0, -60)),
+ ((4, 10), (0, 20)),
+ ((31, 10), (0, 0 )),
+ ((4, 40), (20, 10)),
+ ((5, 40), (20, -15)),
+ ((7, 40), (20, -36)),
+ ((8, 96), (40, -45)),
+ ((9, 96), (40, 0)),
+ ((10, 96), (40, -70)),
+ ((11, 96), (40, -25)),
+ ((12, 96), (40, 30)),
+ ((13, 120), (60, -60)),
+ ((14, 120), (60, -105)),
+ ((15, 120), (60, -85)),
+ ((16, 120), (60, -15)),
+ ((17, 120), (60, 30)),
+ ((18, 120), (60, 5)),
+ ((19, 120), (60, -35)),
+ ((20, 120), (60, 50)),
+ ((14, 132), (80, -160)),
+ ((16, 132), (80, -25)),
+ ((18, 132), (80, 45)),
+ ((22, 132), (80, -55)),
+ ((23, 132), (80, -130)),
+ ((24, 132), (80, -110)),
+ ((25, 132), (80, -100)),
+ ((26, 132), (80, 0)),
+ ((27, 132), (80, 20 )),
+ ((29, 132), (80, -75)),
+ ((30, 132), (80, 80))]
+
+
+ possD = dict([(tp, np.array(pos)) for tp, pos in poss])
+
+ return dict([(g, possD[(g.pid, g.t)]) for g in G.nodes()])
+
+def avgOr(ns, df=0.0):
+ if not ns: return df
+
+ return np.mean(ns)
+
+def avg2(v1, v2):
+ return np.mean([v1, v2])
+
+def ganiso(n):
+ return n.ga
+
+def gr(n):
+ return n.gr
+
+def drawTGraph(G, attr):
+ pos = nodePoss(G)
+ grates = [attr(g) for g in G.nodes()]
+ ws = nx.get_edge_attributes(G, 'connTyp')
+ nodes = nx.draw_networkx_nodes(G, pos, node_color=grates,
+ cmap=plt.cm.Blues)#, vmin=0.0, vmax=0.3)
+
+ labels = {n:str(n.pid) for n in G.nodes()}
+ nx.draw_networkx_labels(G, pos, labels, font_size=10,
+ font_family='sans-serif')#, font_color='slategrey')
+
+ blueEdges = [e for e in G.edges() if ws[e] == 1]
+ redEdges = [e for e in G.edges() if ws[e] == 2]
+ dotEdges = [e for e in G.edges() if ws[e] == 3]
+
+ print(len(blueEdges))
+ print(len(redEdges))
+ print(len(dotEdges))
+
+ nx.draw_networkx_edges(G, pos, edgelist=blueEdges, edge_color='b')
+ nx.draw_networkx_edges(G, pos, edgelist=redEdges, edge_color='r')
+ nx.draw_networkx_edges(G, pos, edgelist=dotEdges, style='dashed')
+
+ plt.colorbar(nodes, shrink=0.5)#, ticks=[0.0, 0.15, 0.3])
+ plt.axis('off')
+
+ plt.savefig("transitionGrowth.svg", dpi=300, box_inches='tight')
diff --git a/growthAnalysis/statesN.py b/growthAnalysis/statesN.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4281845b49ff5ac610e5b77af0021d4296d09ec
--- /dev/null
+++ b/growthAnalysis/statesN.py
@@ -0,0 +1,62 @@
+
+
+def fset(it): return frozenset(it)
+
+states = {fset(['ANT', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY', 'MP']) : 1,
+ fset(['ANT', 'ATML1', 'ETTIN', 'LFY', 'MP', 'PHB_PHV', 'PUCHI',
+ 'REV']) : 2,
+ fset(['ANT', 'ATML1', 'LFY', 'MP', 'PHB_PHV', 'PUCHI', 'REV']) : 3,
+ fset(['ATML1', 'CUC1_2_3', 'MP', 'STM']) : 4,
+ fset(['ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY', 'MP',
+ 'SVP']):5,
+ fset(['ANT', 'AP1', 'AP2', 'LFY', 'MP', 'PHB_PHV', 'PUCHI', 'REV',
+ 'STM', 'SVP']): 6,
+ fset(['ANT', 'AP1', 'AP2', 'ATML1', 'LFY', 'MP', 'PHB_PHV', 'REV',
+ 'STM', 'SVP']): 7,
+ fset(['AHP6', 'ANT', 'AP1', 'AP2', 'ATML1', 'ETTIN', 'LFY', 'MP',
+ 'SEP1', 'SEP2']): 8,
+ fset(['ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY', 'MP',
+ 'SEP1', 'SEP2', 'STM', 'SVP']): 9,
+ fset(['ANT', 'AP1', 'AP2', 'ATML1', 'ETTIN', 'LFY', 'MP', 'PHB_PHV',
+ 'REV', 'SEP1', 'SEP2', 'STM']): 10,
+ fset(['ANT', 'AP1', 'AP2', 'ATML1', 'ETTIN', 'LFY', 'MP', 'SEP1',
+ 'SEP2']): 11,
+ fset(['ANT', 'AP1', 'AP2', 'ATML1', 'ETTIN', 'LFY', 'MP', 'SEP1',
+ 'SEP2', 'STM']): 12,
+ fset(['AG', 'AP3', 'ATML1', 'ETTIN', 'MP', 'PHB_PHV', 'PI', 'REV',
+ 'SEP1', 'SEP2', 'SEP3', 'STM']): 13,
+ fset(['AG', 'ATML1', 'CLV3', 'ETTIN', 'MP', 'PHB_PHV', 'REV', 'SEP1',
+ 'SEP2', 'SEP3', 'STM']): 14,
+ fset(['AG', 'ATML1', 'ETTIN', 'MP', 'PHB_PHV', 'REV', 'SEP1', 'SEP2',
+ 'SEP3', 'STM']): 15,
+ fset(['AHP6', 'ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY',
+ 'MP', 'SEP1', 'SEP2']): 16,
+ fset(['ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY', 'MP',
+ 'SEP1', 'SEP2']): 17,
+ fset(['ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'LFY', 'MP', 'PHB_PHV', 'REV',
+ 'SEP1', 'SEP2']): 18,
+ fset(['AP1', 'AP2', 'AP3', 'ATML1', 'CUC1_2_3', 'MP', 'PHB_PHV', 'PI',
+ 'REV', 'SEP1', 'SEP2', 'STM']): 19,
+ fset(['AP1', 'AP2', 'ATML1', 'CUC1_2_3', 'MP', 'SEP1', 'SEP2', 'STM']): 20,
+ fset(['AP1', 'AP2', 'MP', 'SEP1', 'SEP2', 'STM']): 21,
+ fset(['AG', 'ANT', 'AP3', 'ATML1', 'ETTIN', 'MP', 'PHB_PHV', 'PI', 'SEP1',
+ 'SEP2', 'SEP3', 'STM']): 22,
+ fset(['AG', 'ANT', 'ATML1', 'ETTIN', 'MP', 'PHB_PHV', 'REV', 'SEP1', 'SEP2',
+ 'SEP3', 'STM']): 23,
+ fset(['AG', 'ANT', 'ATML1', 'ETTIN', 'MP', 'PHB_PHV', 'REV', 'SEP1', 'SEP2',
+ 'SEP3', 'STM', 'SUP']): 24,
+ fset(['AG', 'AP3', 'ATML1', 'CUC1_2_3', 'MP', 'PI', 'SEP1', 'SEP2', 'SEP3',
+ 'STM']): 25,
+ fset(['AHP6', 'ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'LFY', 'MP', 'PHB_PHV',
+ 'REV', 'SEP1', 'SEP2']): 26,
+ fset(['ANT', 'AP1', 'AP2', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY',
+ 'SEP1', 'SEP2']): 27,
+ fset(['ANT', 'AP1', 'AP2', 'AS1', 'FIL', 'LFY', 'SEP1', 'SEP2']): 28,
+ fset(['AP1', 'AP2', 'AP3', 'ATML1', 'CUC1_2_3', 'MP', 'PI', 'SEP1',
+ 'SEP2', 'STM']): 29,
+ fset(['AP1', 'AP2', 'AP3', 'ATML1', 'CUC1_2_3', 'MP', 'SEP1', 'SEP2',
+ 'STM']): 30,
+ fset(['ANT', 'AS1', 'ATML1', 'ETTIN', 'FIL', 'LFY', 'MP', 'PUCHI']): 31}
+
+def get_state(gs):
+ return states.get(gs, -1)