Python ConcreteModel.Q方法代码示例

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import Q [as 别名]

#.........这里部分代码省略.........
                    if s >= thr:
                        alleles_I.setdefault(a.name, set()).add(seq)
                    imm[seq, a.name] = min(1., max(0.0, 1.0 - math.log(s, 50000)))
                else:
                    if s > self.__thresh.get(a.name, -float("inf")):
                        alleles_I.setdefault(a.name, set()).add(seq)
                    imm[seq, a.name] = s

            prots = set(pr for pr in p.get_all_proteins())
            cons[seq] = len(prots)
            for prot in prots:
                variations.append(prot.gene_id)
                epi_var.setdefault(prot.gene_id, set()).add(seq)
                var_epi.setdefault(str(seq), set()).add(prot.gene_id)
        self.__peptideSet = peps

        # calculate conservation
        variations = set(variations)
        total = len(variations)
        for e, v in cons.iteritems():
            try:
                cons[e] = v / total
            except ZeroDivisionError:
                cons[e] = 1
        model = ConcreteModel()

        ######################################
        #
        # MODEL DEFINITIONS
        #
        ######################################

        # set definition
        model.Q = Set(initialize=variations)
        model.E = Set(initialize=set(peps.keys()))
        model.TAA = Set(initialize=set(taa))
        model.A = Set(initialize=alleles_I.keys())
        model.G = Set(model.E, initialize=lambda model, e: var_epi[e])
        model.E_var = Set(model.Q, initialize=lambda mode, v: epi_var[v])
        model.A_I = Set(model.A, initialize=lambda model, a: alleles_I[a])

        if self.__included is not None:
            if len(self.__included) > k:
                raise ValueError("More epitopes to include than epitopes to select! "
                                 "Either raise k or reduce epitopes to include.")
        model.Include = Set(within=model.E, initialize=self.__included)

        if overlap > 0:
            def longest_common_substring(model):
                result = []
                for s1,s2 in itr.combinations(model.E,2):
                    if s1 != s2:
                        if s1 in s2 or s2 in s1:
                            result.append((s1,s2))
                        m = [[0] * (1 + len(s2)) for i in xrange(1 + len(s1))]
                        longest, x_longest = 0, 0
                        for x in xrange(1, 1 + len(s1)):
                            for y in xrange(1, 1 + len(s2)):
                                if s1[x - 1] == s2[y - 1]:
                                    m[x][y] = m[x - 1][y - 1] + 1
                                    if m[x][y] > longest:
                                        longest = m[x][y]
                                        x_longest = x
                                else:
                                    m[x][y] = 0
                        if len(s1[x_longest - longest: x_longest]) >= overlap:

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import Q [as 别名]

#.........这里部分代码省略.........
        res_df = results.xs(results.index.values[0][1], level="Method")
        res_df = res_df[res_df.apply(lambda x: any(x[a] > self.__thresh.get(a.name, -float("inf"))
                                                   for a in res_df.columns), axis=1)]

        for tup in res_df.itertuples():
            p = tup[0]
            seq = str(p)
            peps[seq] = p
            for a, s in itr.izip(res_df.columns, tup[1:]):
                if method in ["smm", "smmpmbec", "arb", "comblibsidney"]:
                    try:
                        thr = min(1., max(0.0, 1.0 - math.log(self.__thresh.get(a.name),
                                                      50000))) if a.name in self.__thresh else -float("inf")
                    except:
                        thr = 0

                    if s >= thr:
                        alleles_I.setdefault(a.name, set()).add(seq)
                    imm[seq, a.name] = min(1., max(0.0, 1.0 - math.log(s, 50000)))
                else:
                    if s > self.__thresh.get(a.name, -float("inf")):
                        alleles_I.setdefault(a.name, set()).add(seq)
                    imm[seq, a.name] = s

            prots = set(pr for pr in p.get_all_proteins())
            cons[seq] = len(prots)
            for prot in prots:
                variations.append(prot.gene_id)
                epi_var.setdefault(prot.gene_id, set()).add(seq)
        self.__peptideSet = peps

        #calculate conservation
        variations = set(variations)
        total = len(variations)
        for e, v in cons.iteritems():
            try:
                cons[e] = v / total
            except ZeroDivisionError:
                cons[e] = 1
        model = ConcreteModel()

        #set definition
        model.Q = Set(initialize=variations)

        model.E = Set(initialize=set(peps.keys()))

        model.A = Set(initialize=alleles_I.keys())
        model.E_var = Set(model.Q, initialize=lambda mode, v: epi_var[v])
        model.A_I = Set(model.A, initialize=lambda model, a: alleles_I[a])


        #parameter definition
        model.k = Param(initialize=self.__k, within=PositiveIntegers, mutable=True)
        model.p = Param(model.A, initialize=lambda model, a: probs[a])

        model.c = Param(model.E, initialize=lambda model, e: cons[e],mutable=True)

        #threshold parameters
        model.i = Param(model.E, model.A, initialize=lambda model, e, a: imm[e, a])
        model.t_allele = Param(initialize=0, within=NonNegativeIntegers, mutable=True)
        model.t_var = Param(initialize=0, within=NonNegativeIntegers, mutable=True)
        model.t_c = Param(initialize=0.0, within=NonNegativeReals, mutable=True)

        # Variable Definition
        model.x = Var(model.E, within=Binary)
        model.y = Var(model.A, within=Binary)
        model.z = Var(model.Q, within=Binary)

        # Objective definition
        model.Obj = Objective(
            rule=lambda model: sum(model.x[e] * sum(model.p[a] * model.i[e, a] for a in model.A) for e in model.E),
            sense=maximize)


        #Obligatory Constraint (number of selected epitopes)
        model.NofSelectedEpitopesCov = Constraint(rule=lambda model: sum(model.x[e] for e in model.E) <= model.k)

        #optional constraints (in basic model they are disabled)
        model.IsAlleleCovConst = Constraint(model.A,
                                            rule=lambda model, a: sum(model.x[e] for e in model.A_I[a]) >= model.y[a])
        model.MinAlleleCovConst = Constraint(rule=lambda model: sum(model.y[a] for a in model.A) >= model.t_allele)

        model.IsAntigenCovConst = Constraint(model.Q,
                                             rule=lambda model, q: sum(model.x[e] for e in model.E_var[q]) >= model.z[q])
        model.MinAntigenCovConst = Constraint(rule=lambda model: sum(model.z[q] for q in model.Q) >= model.t_var)
        model.EpitopeConsConst = Constraint(model.E,
                                            rule=lambda model, e: (1 - model.c[e]) * model.x[e] <= 1 - model.t_c)

        #generate instance
        self.instance = model
        if self.__verbosity > 0:
            print "MODEL INSTANCE"
            self.instance.pprint()

        #constraints
        self.instance.IsAlleleCovConst.deactivate()
        self.instance.MinAlleleCovConst.deactivate()
        self.instance.IsAntigenCovConst.deactivate()
        self.instance.MinAntigenCovConst.deactivate()
        self.instance.EpitopeConsConst.deactivate()

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import Q [as 别名]
def create_model(model_name, nodes, links, type_nodes, type_links, timesteps, params):

    m = ConcreteModel(name=model_name)

    # SETS

    # basic sets
    m.Nodes = Set(initialize=nodes)  # nodes
    m.Links = Set(initialize=links)  # links
    m.TS = Set(initialize=timesteps, ordered=True)  # time steps

    # all nodes directly upstream from a node
    def NodesIn_init(m, node):
        retval = []
        for (i, j) in m.Links:
            if j == node:
                retval.append(i)
        return retval

    m.NodesIn = Set(m.Nodes, initialize=NodesIn_init)

    # all nodes directly downstream from a node
    def NodesOut_init(m, node):
        retval = []
        for (j, k) in m.Links:
            if j == node:
                retval.append(k)
        return retval

    m.NodesOut = Set(m.Nodes, initialize=NodesOut_init)

    # sets (nodes or links) for each template type
    for k, v in type_nodes.items():
        exec("m.{} = Set(within=m.Nodes, initialize={})".format(k.replace(" ", "_"), v))
    for k, v in type_links.items():
        exec("m.{} = Set(within=m.Links, initialize={})".format(k.replace(" ", "_"), v))

    # sets for non-storage nodes
    m.NonReservoir = m.Nodes - m.Reservoir
    m.DemandNodes = m.NonReservoir - m.Junction

    # these are collected to initialize the node-block/link-block sets
    demand_node_blocks = []
    reservoir_blocks = []
    link_blocks = []

    # set - all blocks in each demand or reservoir node, and identify node-blocks
    def NodeBlockLookup_init(m, node):
        if "Priority" in params["node"] and node in params["node"]["Priority"]:
            blocks = params["node"]["Priority"][node].columns
        else:
            blocks = [0]  # every node should have a priority
        node_blocks = [(node, b) for b in blocks]
        if node in m.DemandNodes:
            demand_node_blocks.extend(node_blocks)
        elif node in m.Reservoir:
            reservoir_blocks.extend(node_blocks)
        return blocks

    m.DemandNodeBlockLookup = Set(m.DemandNodes, initialize=NodeBlockLookup_init)
    m.ReservoirBlockLookup = Set(m.Reservoir, initialize=NodeBlockLookup_init)

    # set - all blocks in each link
    def LinkBlockLookup_init(m, i, j):
        if "Priority" in params["link"] and (i, j) in params["node"]["Priority"]:
            blocks = params["link"]["Priority"][(i, j)].columns
        else:
            blocks = [0]  # every link should have a priority
            # return Set.End
        link_blocks.extend([(i, j, b) for b in blocks])
        return blocks

    m.LinkBlockLookup = Set(m.Links, initialize=LinkBlockLookup_init)

    # create node-block and link-block sets
    m.DemandNodeBlocks = Set(initialize=demand_node_blocks)
    m.ReservoirBlocks = Set(initialize=reservoir_blocks)
    m.LinkBlocks = Set(initialize=link_blocks)

    # VARIABLES

    m.D = Var(m.DemandNodes * m.TS, domain=NonNegativeReals)  # delivery to demand nodes
    m.D_DB = Var(m.DemandNodeBlocks * m.TS, domain=NonNegativeReals)  # delivery to demand nodes
    m.D_surplus = Var(m.DemandNodes * m.TS, domain=NonNegativeReals)  # delivery to demand nodes
    m.S = Var(m.Reservoir * m.TS, domain=NonNegativeReals)  # storage
    m.S_RB = Var(m.ReservoirBlocks * m.TS, domain=NonNegativeReals)  # storage
    m.S_surplus = Var(m.Reservoir * m.TS, domain=NonNegativeReals)  # storage

    m.G = Var(m.Nodes * m.TS, domain=NonNegativeReals)  # gain (local inflow)
    m.L = Var(m.Nodes * m.TS, domain=NonNegativeReals)  # loss (local outflow)
    m.I = Var(m.Nodes * m.TS, domain=NonNegativeReals)  # total inflow to a node
    m.O = Var(m.Nodes * m.TS, domain=NonNegativeReals)  # total outflow from a node

    m.Q = Var(m.Links * m.TS, domain=NonNegativeReals)  # flow in links
    m.Q_LB = Var(m.LinkBlocks * m.TS, domain=NonNegativeReals)  # flow in links
    m.Q_surplus = Var(m.Links * m.TS, domain=NonNegativeReals)  # flow in links

    # PARAMETERS

    # IMPORTANT: Defaults should not be set here
#.........这里部分代码省略.........