growth_control.py

import os
import math
from itertools import permutations, product, repeat, takewhile
from functools import partial, reduce
from operator import add, itemgetter
from collections import defaultdict
from subprocess import call
from typing import Dict, List, Tuple, Set

import ast
import astor
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib as mpl
from scipy.signal import argrelmax
from scipy.stats import gaussian_kde
from scipy.stats import ttest_ind
from more_itertools import interleave_longest
from mpl_toolkits.mplot3d import Axes3D

import common.statesN as f
import common.seg as seg
import common.edict as edict
import common.lin as lin
from common.seg import STissue
import grates as grs
import patternTransitions as p

geneNms: List[str] = ['"AG"',
                      '"AHP6"',
                      '"ANT"',
                      '"AP1"',
                      '"AP2"',
                      '"AP3"',
                      '"AS1"',
                      '"ATML1"',
                      '"CLV3"',
                      '"CUC1_2_3"',
                      '"ETTIN"',
                      '"FIL"',
                      '"LFY"',
                      '"MP"',
                      '"PHB_PHV"',
                      '"PI"',
                      '"PUCHI"',
                      '"REV"',
                      '"SEP1"',
                      '"SEP2"',
                      '"SEP3"',
                      '"STM"',
                      '"SUP"',
                      '"SVP"',
                      '"WUS"']

mpl.rcParams.update(mpl.rcParamsDefault)
sts_ts = f.states_per_t()

fst = itemgetter(0)
snd = itemgetter(1)


def set_plot_params(fontsize=11):
    sns.set_style("whitegrid", {'grid.linestyle': '--'})
    plt.rc('xtick', labelsize=fontsize)
    plt.rc('ytick', labelsize=fontsize)
    plt.rc('axes', labelsize=fontsize)
    params = {'legend.fontsize': fontsize}
    plt.rcParams.update(params)


class BExpr2():

    def __init__(self, e_str):
        self.e = ge(e_str)

    def __hash__(self):
        gns = sorted([hash(g) for g in getGenes(getTopExpr(self.e))])
        op = hash(ast.dump(getTopExpr(self.e).op))

        return sum(gns) + op

    def __eq__(self, other):
        return hash(self) == hash(other)

    def __repr__(self):
        import astor
        return astor.to_source(self.e).replace('"', '').strip()


class BExpr2Region():

    def __init__(self, t, ts, e_str):
        self.t = t
        self.e = ge(e_str)
        self.ts = ts

    def __hash__(self):
        return sum(get_sts_expr(self.ts, self.e, sts_ts[self.t]))

    def __eq__(self, other):
        return hash(self) == hash(other)


class GPairRegions():
    def __init__(self, t, ts, g1, g2):
        self.t = t
        self.g1 = g1
        self.g2 = g2
        self.e1 = ge("'{g1}' and '{g2}'".format(g1=g1, g2=g2))
        self.e2 = ge("'{g1}' and not '{g2}'".format(g1=g1, g2=g2))
        self.ts = ts

    def __hash__(self):
        sts = sts_ts[self.t]
        sts_e1 = get_sts_expr(self.ts, self.e1, sts)
        sts_e2 = get_sts_expr(self.ts, self.e2, sts)

        return (hash(frozenset(sts_e1)) +
                hash(frozenset(sts_e2)))

    def __eq__(self, other):
        return hash(self) == hash(other)

    def plot(self, txt=""):
        d = dict()
        grs_state = get_grs_state()
        sts = sts_ts[self.t]
        sts_e1 = get_sts_expr(self.ts, self.e1, sts)
        sts_e2 = get_sts_expr(self.ts, self.e2, sts)
        gr_e1 = meanor0(reduce(add, [grs_state[st] for st in sts_e1], []))
        gr_e2 = meanor0(reduce(add, [grs_state[st] for st in sts_e2], []))

        he, le = [e for e, gr in sorted([(self.e1, gr_e1), (self.e2, gr_e2)],
                                        key=snd, reverse=True)]

        for c in self.ts:
            if seg.evalB(getTopExpr(he), c.exprs):
                d[c.cid] = 1.0
            else:
                d[c.cid] = 0.5

        for c in self.ts:
            if seg.evalB(getTopExpr(le), c.exprs):
                d[c.cid] = 0.0

        f = plot_ts_q_(self.ts, d,
                        lb="{g1}-{g2}".format(g1=self.g1, g2=self.g2),
                        bounds=(0.1, 0.9), txt=txt)

        return f, d

    def get_regions(self) -> Tuple[Set[int], Set[int]]:
        grs_state = get_grs_state()
        sts = sts_ts[self.t]
        sts_e1 = get_sts_expr(self.ts, self.e1, sts)
        sts_e2 = get_sts_expr(self.ts, self.e2, sts)
        gr_e1 = meanor0(reduce(add, [grs_state[st] for st in sts_e1], []))
        gr_e2 = meanor0(reduce(add, [grs_state[st] for st in sts_e2], []))

        he, le = [e for e, gr in sorted([(self.e1, gr_e1), (self.e2, gr_e2)],
                                        key=snd, reverse=True)]

        r1 = list()
        r2 = list()
        for c in self.ts:
            if seg.evalB(getTopExpr(he), c.exprs):
                r1.append(c.cid)
            elif seg.evalB(getTopExpr(le), c.exprs):
                r2.append(c.cid)

        return set(r1), set(r2)

    def plot_distr(self):
        set_plot_params(fontsize=12)
        sns.set_style("white")
        cls = sns.color_palette("coolwarm")
        grs_state = get_grs_state()
        sts = sts_ts[self.t]
        sts_e1 = get_sts_expr(self.ts, self.e1, sts)
        sts_e2 = get_sts_expr(self.ts, self.e2, sts)
        grs_e1 = reduce(add, [grs_state[st] for st in sts_e1])
        grs_e2 = reduce(add, [grs_state[st] for st in sts_e2])

        gr_e1 = np.mean(reduce(add, [grs_state[st] for st in sts_e1]))
        gr_e2 = np.mean(reduce(add, [grs_state[st] for st in sts_e2]))

        h_grs, l_grs = [e for e, gr in sorted([(grs_e1, gr_e1), (grs_e2, gr_e2)],
                                              key=snd, reverse=True)]

        he, le = [e for e, gr in sorted([(self.e1, gr_e1), (self.e2, gr_e2)],
                                        key=snd, reverse=True)]

        fig = plt.figure(figsize=(3.5, 3.7))
        ax = fig.add_subplot('111')

        ax.set_xlim(-0.05, 0.18)
        
        sns.distplot(h_grs, bins=10, color=cls[-1], kde=True, ax=ax,
                     label=astor.to_source(he).replace("'", ""))
        sns.distplot(l_grs, bins=10, color=cls[0], kde=True, ax=ax,
                     label=astor.to_source(le).replace("'", ""))
        ax.set_xlabel(r'$\mu$m/h')
        #plt.legend(fontsize=10, loc='upper right', frameon=False)
        gns = "{g1}-{g2}".format(g1=self.g1, g2=self.g2)
        print(gns, gratio(h_grs, l_grs))
        plt.savefig("{g1}-{g2}-regions.png".format(g1=self.g1, g2=self.g2), dpi=300)
        plt.show()

        return

    def toPair(self):
        return (self.g1, self.g2)


def groupByRegion(t, ts, es: List[str]) -> Dict[int, List[str]]:
    res = defaultdict(list)

    regs = set([hash(BExpr2Region(t, ts, e)) for e in es])
    reg_ids = {reg:i for i, reg in enumerate(regs)}

    for e in es:
        e_rid = reg_ids[hash(BExpr2Region(t, ts, e))]
        res[e_rid].append(e)

    return res


def mean_gratio(ps):
    return meanor0([abs(v) for p, v in ps])


def mean_gratio_(ps):
    return meanor0([v for p, v in ps])


def groupByRegionGPairs(t: int,
                        ts: STissue,
                        ress: List[Tuple[Tuple[str, str], float]],
                        n=6,
                        prc=0.5):
    d = defaultdict(list)
    ps_nnan = [(GPairRegions(t, ts, p[0], p[1]), v) for p, v in ress
               if not math.isnan(v)]
    nps = int(np.ceil(len(ps_nnan)*prc))
    ps = sorted(ps_nnan,
                key=lambda x: abs(x[1]),
                reverse=True)[:nps]
    regs = set([hash(p) for p, v in ps])

    for p, v in ps:
        d[hash(p)].append((p.toPair(), v))

    d_ = sorted([(k, v) for k, v in d.items()],
                key=lambda ps: mean_gratio(ps[1]),
                reverse=True)[:n]

    return [(i+1, v) for i, (k, v) in zip(range(len(d_)), d_)]


def plotRegions(t, ts, d):
    fns = list()
    for rid, ps in d:
        r = GPairRegions(t, ts, ps[0][0][0], ps[0][0][1])
        f = r.plot(txt="{rid} ({gr:.3f})".format(rid=str(rid),
                                                 gr=mean_gratio(ps)))
        fns.append(f)

    montage_fimgs(fns, im_label="t{t}_regions".format(t=t))

    for fn in fns:
        os.remove(fn)

    return "t{t}_regions.png".format(t=t)


def cluster_states(grs_state: Dict[int, List[float]]):
    mgrs = {st:np.mean(grs) for st, grs in grs_state.items()}
    g = gaussian_kde(np.array(list(mgrs.values())))

    xs = np.arange(-0.05, 0.2, 0.001)
    ys = g.evaluate(xs)

    cms = argrelmax(ys)[0]

    xsm = np.array(xs[cms])
    ysm = np.array(ys[cms])

    cls = {st:np.argmin(np.abs(xsm - grm))
           for st, grm in mgrs.items()}

    return g, xsm, ysm, cls


def plot_clustering(g, xsm, ysm):
    xs = np.arange(-0.05, 0.2, 0.001)
    ys = g.evaluate(xs)
    cls = sns.color_palette()

    plt.plot(xs, ys, color=cls[0])
    plt.plot(xsm, ysm, 'o', color=cls[1])
    plt.xlabel("growth rate")
    plt.ylabel("density estimate")
    plt.savefig("clusters_grates_states.png", dpi=300)
    plt.show()


def get_grs_state() -> Dict[int, List[float]]:
    tss, linss = lin.mkSeries1(d="../data/FM1/tv/",
      	                       dExprs="../data/geneExpression/",
			       linDataLoc="../data/FM1/tracking_data/",
			       ft=lambda t: t in {10, 40, 96, 120, 132})
    lin.filterL1_st(tss)
    G = p.mkTGraphN()

    grates = grs.grates_avg_cons(tss, linss)
    grates_pats = p.addCombPatterns(grates, tss)
    grates_state = p.getGAnisosPerPattern(G, grates_pats)

    return grates_state


def ge(s: str): return ast.parse(s)


def construct_expr(nms, fns):
    return ge(" ".join(interleave_longest(nms, fns)))


def mk_bexprs(k, geneNms):
    not_genes = ["not {gn}".format(gn=gn) for gn in geneNms]
    not_gene_nms = geneNms + not_genes

    fnss = [['and', 'or']] * (k-1)

    for gns in permutations(not_gene_nms, k):
        for fns in product(*fnss):
            bexpr = construct_expr(gns, fns)
            yield bexpr

    return


def search_bexprs(ts, k, obj):
    ress = dict()
    for i, bexpr in enumerate(mk_bexprs(k, geneNms)):
        # cids = ts.filterGBExpr(bexpr)
        ress[bexpr] = obj(cids2=bexpr)

    return ress


def hash_expr(ts, e):
    return np.sum(ts.filterGBExpr(e))


def get_cids_highg(ts):
    sts = f.statesI()
    grs_state = get_grs_state()
    g, xsm, ysm, cls = cluster_states(grs_state)
    stsHigh = [st for st, cl in cls.items() if cl == 1]

    cids = reduce(add, [ts.filterGs(sts[st]) for st in stsHigh])

    return cids


def get_cids_sts(ts, sts):
    st_genes = f.statesI()
    return reduce(add, [ts.filterGs(st_genes[st]) for st in sts])


def bacc(ts, cids1, cids2):
    scids1 = set(cids1)
    scids2 = set(cids2)

    pos = sum([1 for c in ts if c.cid in scids1])
    neg = len(list(ts)) - pos

    tp = sum([1 for c in ts
              if c.cid in scids1 and c.cid in scids2])
    tn = sum([1 for c in ts
              if c.cid not in scids1 and c.cid not in scids2])

    if not pos == 0:
        tpr = tp / pos
    else:
        tpr = 1.0

    if not neg == 0:
        tnr = tn / neg
    else:
        tnr = 1.0

    return 0.5*(tpr + tnr)


def meanor0(xs):
    if xs:
        return np.mean(xs)
    else:
        return 0.0


def avg_grate(grates, cids2):
    grs_cids = [grates.get(cid, None) for cid in cids2]
    grs_cids_ = [gr for gr in grs_cids if gr]

    return meanor0(grs_cids_)


def gratio(xs, xs1):
    if not xs or not xs1:
        return 0.0

    m1 = np.mean(xs)
    m2 = np.mean(xs1)

    return (m1 - m2) / (m1 + m2)


def states_d(ts, cids2, sts, grs_state_m):
    sts_envs = f.statesI_env(seg.geneNms)
    e = getTopExpr(cids2)
    mgrs_sts_t = [grs_state_m[st] for st in sts
                  if seg.evalB(e, sts_envs[st])]
    mgrs_sts_f = [grs_state_m[st] for st in sts
                  if not seg.evalB(e, sts_envs[st])]

    return gratio(mgrs_sts_t, mgrs_sts_f)


def go_search(t, k):
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)

    ts = tss[t]
    cids_pat = get_cids_highg(ts)
    bacc_obj = partial(bacc, ts=ts, cids1=cids_pat)

    ress = search_bexprs(ts, k, bacc_obj)

    return ress

def go_search1(t, k):
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)
    ts = tss[t]

    grates = grs.grates_avg_cons(tss, linss)
    gr_obj = partial(avg_grate, grates=grates[t])

    ress = search_bexprs(ts, k, gr_obj)

    return ress

def get_scores(t, es_reg, reg_sts):
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)
    ts = tss[t]
    sts = f.states_per_t()[t]

    scores = defaultdict(list)
    for ri, ess in es_reg.items():
        for es1 in ess:
            r_sts = reg_sts[ri]
            cids1 = get_cids_sts(ts, r_sts)
            cids2 = ts.filterGBExpr(ge(es1))
            s = bacc(ts, cids1, cids2)
            scores[ri].append((es1, s))

    return scores

def filterNonEmpty(ts, es):
    #an expression makes sense
    #if all the genes are expressed in the tissue
    es_ = list()
    for e in es:
        gs = ["'{g}'".format(g=g.replace("not", "").strip())
              for g in getGenes(getTopExpr(ge(e)))]
        gs_on = [g for g in gs if ts.filterGBExpr(ge(g))]
        if len(gs) == len(gs_on): es_.append(e)

    return es_

def const(cids2, y):
    return y

def go_search2(t, k):
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)
    ts = tss[t]
    sts = f.states_per_t()[t]

    grs_state = get_grs_state()
    mgrs_state = {st:np.mean(grs)
                  for st, grs in grs_state.items()}
    st_gr_obj = partial(states_d,
                        ts=ts,
                        sts=sts,
                        grs_state_m=mgrs_state)

    ress = search_bexprs(ts, k, st_gr_obj)
    ress_ = pprint_list(ress)

    max_v = ress_[0][1]
    min_v = ress_[-1][1]
    ress_max_elems = list(takewhile(lambda p: p[1] == max_v, ress_))
    ress_min_elems = list(takewhile(lambda p: p[1] == min_v, reversed(ress_)))

    es = [e for e, v in ress_max_elems + ress_min_elems]
    es_reg = groupByRegion(t, ts, es)

    for rid, es in es_reg.items():
        es_ = [repr(e) for e in set([BExpr2(e) for e in es])]
        with open("reg{i}_t{t}_exprs.txt".format(i=rid, t=t), "w+") as fout:
            fout.write("\n".join(es_))

    plot_state_mgrs(t)
    reg_sts = {rid:get_sts_expr(ts, ge(es[0]), sts)
               for rid, es in es_reg.items()}

    for rid, es in es_reg.items():
        plot_ts_binary(ts, ts.filterGBExpr(ge(es[0])),
                       "reg{i}_t{t}".format(i=rid, t=t),
                       txt="states: " + ", ".join(list(map(str, reg_sts[rid]))))

    es_reg_ = dict()
    for rid, es in es_reg.items():
        es_reg_[rid] = filterNonEmpty(ts, list(repr(e) for e in set([BExpr2(e) for e in es])))


    return ress, es_reg_, reg_sts

def go_search3(t, k):
    tss, linss = lin.mkSeries()
    lin.filterL1(tss)
    ts = tss[t]
    sts = f.states_per_t()[t]

    def_obj = partial(const, y=0.0)
    ress = search_bexprs(ts, k, def_obj)
    ress_ = pprint_list(ress)

    es = [e for e, v in ress_]
    es_reg = groupByRegion(t, ts, es)

    reg_sts = {rid:get_sts_expr(ts, ge(es[0]), sts)
               for rid, es in es_reg.items()}

    return es_reg, reg_sts

def invert_dlist(d):
    d_ = dict()
    for k, xs in d.items():
        for x in xs:
            d_[x] = k

    return d_

def vis_stripplot(d):
    fig = plt.figure()
    ax = fig.add_subplot('111')
    sns.set_style("whitegrid", {'grid.linestyle':'--'})
    xs_labs = sorted(d.keys())
    ax.yaxis.set_major_locator(ticker.MultipleLocator(0.05))

    xs = reduce(add, [list(repeat(x, len(d[x])))
                      for i, x in enumerate(xs_labs)])
    ys = reduce(add, [d[x] for x in xs_labs])
    cls = sns.color_palette()
    sns.stripplot(xs, ys, palette=[cls[0]], size=5, alpha=0.75,
                  jitter=False)
    ax.set_xticks(list(range(len(xs_labs))))
    ax.set_ylabel("BAcc")
    ax.set_xlabel("expr length")

    return ax

def vis_stripplot_hue2(d, i):
    fig = plt.figure(figsize=(12, 6))
    ax = fig.add_subplot('111')
    sns.set_style("whitegrid", {'grid.linestyle':'--'})

    xs_labs = sorted(list(set([gn.replace("not", "").strip()
                               for gn in d.keys()])))

    xs = reduce(add, [list(repeat(x, len(d[x])*2))
                      for i, x in enumerate(xs_labs)])
    hs = reduce(add, [(list(repeat("g", len(d[x]))) +
                       list(repeat("not g", len(d[x]))))
                      for i, x in enumerate(xs_labs)])
    ys = reduce(add, [d[x] + d["not {x}".format(x=x)]
                      for i, x in enumerate(xs_labs)])

    cls = sns.color_palette('Blues')
    cls_not = sns.color_palette('Reds')
    sns.stripplot(x=xs, y=ys, hue=hs, palette=[cls[i], cls_not[i]],
                  size=4, alpha=0.6, dodge=True)
    ax.set_xticks(list(range(len(xs_labs))))
    ax.set_ylim((0.01, 0.11))
    plt.xticks(rotation=90)

    plt.savefig("grates_gene.png", dpi=300)

    plt.show()


def cellsToDat(fn, cells):
    # so we can use newman
    nvars = len(cells[0].exprs.keys()) + 5
    ncells = len(cells)
    ntpoints = 1

    header = "\n".join([str(ntpoints),
                        " ".join([str(ncells), str(nvars), "0"])])

    with open(fn, 'w+') as fout:
        fout.write("\n".join([header] + [cell.toOrganism() for cell in cells]))


def plot_ts_binary(ts_, cids, lb="region", txt=""):
    from copy import deepcopy
    ts = deepcopy(ts_)

    def toNewmanFn(fn):
        (fn, ext) = os.path.splitext(fn)
        return fn + "000" + ".tif"

    def toNewmanFnPng(fn):
        (fn, ext) = os.path.splitext(fn)
        return fn + "000" + ".png"

    def ind(ts, g):
        geneNms = list(ts)[0].geneNms
        return geneNms.index(g)

    convertCmd = "convert"
    confFile = "/Users/s1437009/Organism/tools/plot/bin/newmanInit.conf"
    visCmd = "/Users/s1437009/Organism/tools/plot/bin/newman"

    lbGene = "ANT"
    scids = set(cids)
    for c in ts:
        if c.cid in scids:
            c.exprs[lbGene] = True
        else:
            c.exprs[lbGene] = False

    cellsToDat(lb + ".data", list(ts))

    call([visCmd,
          "-shape", "sphereVolume",
          "-d", "3",
          lb + ".data",
          "-column", str(ind(ts, lbGene) + 1),
          "-schema", str(0),
          "-output", "tiff",
          lb,
          "-camera", confFile,
          "-size", str(720),
          "-min", str(-0.03),
          "-max", str(1.03)])

    call([convertCmd,
          toNewmanFn(lb),
          "-fill", "white",
          "-font", "Times-New-Roman",
          "-pointsize", str(30),
          "-undercolor", "'#85929E'",
          "-gravity", "North",
          "-annotate", "+0+5",
          "' {txt} '".format(txt=txt.replace("'", "")),
          "{lb}.png".format(lb=lb)])

    os.remove(lb+".data")
    os.remove(toNewmanFn(lb))

    call(["open",
          "{lb}.png".format(lb=lb)])


def plot_ts_q_(ts, d, lb="vals", bounds=(0, 1), txt=""):
    vmin, vmax = bounds

    xs = [c.pos.x for c in ts]
    ys = [c.pos.y for c in ts]
    zs = [c.pos.z for c in ts]
    cs = [d[c.cid] for c in ts]

    fig = plt.figure(figsize=(1.5, 1.5))
    ax = fig.add_subplot(111, projection='3d')
    ax.view_init(elev=104, azim=-89)
    ax.scatter(xs, ys, zs, c=cs, cmap="coolwarm", alpha=0.9, s=7, vmin=vmin,
               vmax=vmax, edgecolors='black', linewidths=0.5)
    ax.set_title(txt)
    plt.axis('off')

    fout = "{lb}.png".format(lb=lb)
    plt.savefig(fout, dpi=300, bbox_inches='tight')

    plt.show()


def plot_distrs(grs, r1, r2, lb):
    grs1 = [grs.get(cid, np.nan) for cid in r1]
    grs2 = [grs.get(cid, np.nan) for cid in r2]

    grs1_ = [gr for gr in grs1 if not math.isnan(gr)]
    grs2_ = [gr for gr in grs2 if not math.isnan(gr)]

    cls = sns.color_palette("coolwarm")

    fig = plt.figure(figsize=(3.5, 3.7))
    ax = fig.add_subplot('111')

    ax.set_xlim(-0.05, 0.18)

    sns.distplot(grs1_, bins=10, color=cls[-1], kde=True, ax=ax,
                 label="region1")
    sns.distplot(grs2_, bins=10, color=cls[0], kde=True, ax=ax,
                 label="region2")
    ax.set_xlabel(r'$\mu$m/h'.format())
    plt.savefig("zone{i}.png".format(i=lb), dpi=300, bbox_inches='tight')
    plt.show()


def plot_ts_q(ts_, d, lb="vals", bounds=(0, 1), txt="",
              interactive=False, size=100, cm=1):
    from copy import deepcopy
    ts = deepcopy(ts_)

    def toNewmanFn(fn):
        (fn, ext) = os.path.splitext(fn)
        return fn + "000" + ".tif"

    def ind(ts, g):
        geneNms = list(ts)[0].geneNms
        return geneNms.index(g)

    confFile = "/Users/s1437009/Organism/tools/plot/bin/newmanInit.conf"
    visCmd = "/Users/s1437009/Organism/tools/plot/bin/newman"
    convertCmd = "convert"

    lbGene = "ANT"
    for c in ts:
        if c.cid in d:
            c.exprs[lbGene] = d[c.cid]

    cellsToDat(lb + ".data", list(ts))

    call([visCmd,
          "-shape", "sphereVolume",
          "-d", "3",
          lb + ".data",
          "-column", str(ind(ts, lbGene) + 1),
          "-schema", str(cm),
          "-output", "tiff",
          lb,
          "-min", str(bounds[0]),
          "-max", str(bounds[1]),
          "-size", str(size),
          "-camera", confFile])

    if lb != "":
        call([convertCmd,
              toNewmanFn(lb),
              "-fill", "white",
              "-font", "Times-New-Roman",
              "-pointsize", str(15),
              "-undercolor", "#696969",
              "-gravity", "North",
              "-annotate", "+0+0",
              " {txt} ".format(txt=txt),
              "{lb}.png".format(lb=lb)])
    else:
        call([convertCmd,
              toNewmanFn(lb),
              "{lb}.png".format(lb=lb)])

    if interactive:
        call(["open",
              "{lb}.png".format(lb=lb)])

    os.remove(lb+".data")
    os.remove(toNewmanFn(lb))

    return "{lb}.png".format(lb=lb)

def montage_fimgs(fns, im_label="montage"):
    montageCmd = "montage"
    n = len(fns)
    m = np.ceil(n / 3)
    k = min(n, 3)

    #collage all the images
    call([montageCmd] +
         fns +
         ["-tile",
          str(k) + "x" + str(m),
          "-geometry", "+20+20",
          im_label + ".png"])

def mark_state_gr(ts, grs_state):
    st_genes = f.states
    mgrs = dict([(st, np.mean(grs)) for st, grs in grs_state.items()])
    d = dict()
    for c in ts:
        d[c.cid] = mgrs.get(st_genes.get(c.getOnGenes(), None), 0.0)

    return d

def plot_state_mgrs(t):
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)
    ts = tss[t]

    grs_state = get_grs_state()
    d = mark_state_gr(ts, grs_state)

    plot_ts_q(ts, d, "grates_sts_{t}h".format(t=t))


def get_st_cids_expr(ts, e, sts: List[int]):
    sts_envs = f.statesI_env(seg.geneNms)
    sts_genes = f.statesI()
    e = getTopExpr(e)

    sts_t = [st for st in sts if seg.evalB(e, sts_envs[st])]
    cids_t = set(reduce(add, [ts.filterGs(sts_genes[st]) for st in sts_t]))

    return cids_t

def get_sts_expr(ts, e, sts: List[int]) -> List[int]:
    sts_envs = f.statesI_env(seg.geneNms)
    sts_genes = f.statesI()
    e = getTopExpr(e)

    sts_t = [st for st in sts if seg.evalB(e, sts_envs[st])]

    return sts_t

def showR(r):
    e, gr = r
    e = e.replace("'", "")

    return "{e}, {gr:.2f}".format(e=e, gr=gr)

def getGenesOp(e):
    if type(e.op) is ast.And:
        return reduce(add, [getGenes(e1) for e1 in e.values])
    elif type(e.op) is ast.Or:
        return reduce(add, [getGenes(e1) for e1 in e.values])
    elif type(e.op) is ast.Not:
        return ["not {gn}".format(gn=reduce(add, getGenes(e.operand)))]

def getGenes(e):
    if type(e) is ast.UnaryOp or type(e) is ast.BoolOp:
        return getGenesOp(e)
    elif type(e) is ast.Str:
        return [e.s]
    elif type(e) is ast.Name:
        return [e.id]

def getTopExpr(e):
    return e.body[0].value

def expr(r):
    return r[0]

def val(r):
    return r[1]

def getResPerGene(ress):
    ress_gene = defaultdict(list)

    for expr, res in ress.items():
        for gn in getGenes(getTopExpr(expr)):
            ress_gene[gn].append(res)

    return ress_gene

def pprint_list(ress):
    import astor
    ress_ = sorted([(astor.to_source(k).replace('"', '').strip(), v)
                    for k, v in ress.items()],
                   key=lambda x: x[1],
                   reverse=True)

    return ress_

def getLinGBExpr(tss, e):
    return {t: ts.filterGBExpr(e) for t, ts in tss.items()}

def getValsT(cidsT, valsT, tpoints=set([10, 40, 96, 120, 132])):
    return {t:edict.gets(valsT[t], cids) for t, cids in cidsT.items()
            if t in tpoints}

def gratio(xs, xs1):
    if not xs or not xs1:
        return np.nan

    m1 = np.mean(xs)
    m2 = np.mean(xs1)

    return (m1 - m2) / (m1 + m2)

def fgr(tss, grates, g1, g2):
    bexpr1 = "'{g1}' and not '{g2}'".format(g1=g1, g2=g2)
    bexpr2 = "'{g1}' and '{g2}'".format(g1=g1, g2=g2)

    cidsT = getLinGBExpr(tss, seg.ge(bexpr2))
    cidsT1 = getLinGBExpr(tss, seg.ge(bexpr1))

    grs = getValsT(cidsT, grates)
    grs1 = getValsT(cidsT1, grates)

    return {t:gratio(grs[t], grs1[t])
            for t in grs.keys()}


def get_vals_exprs(tss, e1, e2, vals):
    cidsT = getLinGBExpr(tss, seg.ge(e1))
    cidsT1 = getLinGBExpr(tss, seg.ge(e2))

    grs = getValsT(cidsT, vals)
    grs1 = getValsT(cidsT1, vals)

    return grs, grs1
    
def rgd_gene(tss, g, grates):
    bexpr1 = "'{g1}'".format(g1=g)
    bexpr2 = "not '{g1}'".format(g1=g)

    grs, grs1 = get_vals_exprs(tss, bexpr1, bexpr2, grates)
    
    return {t:gratio(grs[t], grs1[t])
            for t in grs.keys()}

def rgd_gene_(tss, g, grates):
    bexpr1 = "'{g1}'".format(g1=g)
    bexpr2 = "not '{g1}'".format(g1=g)

    grs_t, grs1_t = get_vals_exprs(tss, bexpr1, bexpr2, grates)

    return gratio(reduce(add, grs_t.values()),
                  reduce(add, grs1_t.values()))

def rgd_state(st, grs_state):
    grs_ = reduce(add, [grs for st_id, grs in grs_state.items()
                        if st_id != st])
    gr = grs_state[st]

    return gratio(gr, grs_)

def pval_state(st, grs_state):
    grs_ = reduce(add, [grs for st_id, grs in grs_state.items()
                        if st_id != st])
    gr = grs_state[st]

    _, pval = ttest_ind(grs_, gr)

    return pval

def pvals_all_states(grs_state):
    return {st_id:pval_state(st_id, grs_state)
            for st_id in grs_state.keys()}

def rgd_all_states(grs_state):
    return {st_id:rgd_state(st_id, grs_state)
            for st_id in grs_state.keys()}

def go_single(tss, geneNms):
    d = {10:  dict(),
         40:  dict(),
         96:  dict(),
         120: dict(),
         132: dict()}
    
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)
    grates = grs.grates_avg_cons(tss, linss)

    for g in geneNms:
        v = rgd_gene(tss, g, grates)
        for t, gr in v.items():
            d[t][g] = [gr]

    return d

def go_single_(tss, geneNms):
    d = dict()
    tss, linss = lin.mkSeries()
    lin.filterL1_st(tss)
    grates = grs.grates_avg_cons(tss, linss)

    for g in geneNms:
        v = rgd_gene_(tss, g, grates)
        d[g] = [v]