Python ConcreteModel.p方法代码示例

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
    def test_mutable_novalue_param_equality(self):
        model = ConcreteModel()
        model.x = Var()
        model.p = Param(mutable=True)
        model.p.value = None

        model.c = Constraint(expr=model.x - model.p == 0)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x == model.p)
        self.assertTrue(model.c.upper is model.p)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x + 1 == model.p)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x + 1 == (model.p + 1)**2)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x == model.p + 1)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

        model.c = Constraint(expr=inequality(model.p, model.x,  model.p))
        self.assertTrue(model.c.upper is model.p)
        # GH: Not sure if we are supposed to detect equality
        #     in this situation. I would rather us not, for
        #     the sake of making the code less complicated.
        #     Either way, I am not going to test for it here.
        #self.assertEqual(model.c.equality, <blah>)
        model.del_component(model.c)

        model.c = Constraint(expr=(model.x, model.p))
        self.assertTrue(model.c.upper is model.p)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

        model.c = Constraint(expr=(model.p, model.x))
        self.assertTrue(model.c.upper is model.p)
        self.assertEqual(model.c.equality, True)
        model.del_component(model.c)

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
    def test_compute_time_stage(self):
        model = ConcreteModel()
        model.x = Var()
        model.y = Var([0,1])
        model.z = Var()
        model.p = Param(mutable=True)
        model.cost = Expression([0,1])
        model.cost[0] = model.x + model.y[0]
        model.cost[1] = model.p + model.y[1] + model.y[0]*model.p
        model.o = Objective(expr=model.cost[0] + model.cost[1])
        model.c = ConstraintList()
        model.c.add(model.x >= 1)               # 1
        model.c.add(model.y[0] >= 1)            # 2
        model.c.add(model.p * model.y[0] >= 1)  # 3
        model.c.add(model.y[0] >= model.p)      # 4
        model.c.add(model.p <= model.y[1])      # 5
        model.c.add(model.y[1] <= 1)            # 6
        model.c.add(model.x >= model.p)         # 7
        model.c.add(model.z == 1)               # 8

        model.varstage = VariableStageAnnotation()
        model.varstage.declare(model.x, 1)
        model.varstage.declare(model.y[0], 1, derived=True)
        model.varstage.declare(model.y[1], 2)
        model.varstage.declare(model.z, 2, derived=True)

        model.stagecost = StageCostAnnotation()
        model.stagecost.declare(model.cost[0], 1)
        model.stagecost.declare(model.cost[1], 2)

        model.stochdata = StochasticDataAnnotation()
        model.stochdata.declare(model.p,
                                distribution=UniformDistribution(0, 1))
        sp = EmbeddedSP(model)

        #
        # check variables
        #
        self.assertEqual(sp.compute_time_stage(model.x),
                         min(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.x,
                                               derived_last_stage=True),
                         min(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.y[0]),
                         min(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.y[0],
                                               derived_last_stage=True),
                         max(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.y[1]),
                         max(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.y[1],
                                               derived_last_stage=True),
                         max(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.z),
                         max(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.z,
                                               derived_last_stage=True),
                         max(sp.time_stages))

        #
        # check constraints
        #
        self.assertEqual(sp.compute_time_stage(model.c[1]),
                         min(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.c[1],
                                               derived_last_stage=True),
                         min(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.c[2]),
                         min(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.c[2],
                                               derived_last_stage=True),
                         max(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.c[3]),
                         max(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.c[3],
                                               derived_last_stage=True),
                         max(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.c[4]),
                         max(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.c[4],
                                               derived_last_stage=True),
                         max(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.c[5]),
                         max(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.c[5],
                                               derived_last_stage=True),
                         max(sp.time_stages))

        self.assertEqual(sp.compute_time_stage(model.c[6]),
                         max(sp.time_stages))
        self.assertEqual(sp.compute_time_stage(model.c[6],
                                               derived_last_stage=True),
                         max(sp.time_stages))
#.........这里部分代码省略.........

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
    def test_mutable_novalue_param_upper_bound(self):
        model = ConcreteModel()
        model.x = Var()
        model.p = Param(mutable=True)
        model.p.value = None

        model.c = Constraint(expr=model.x - model.p <= 0)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x <= model.p)
        self.assertTrue(model.c.upper is model.p)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x + 1 <= model.p)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x <= model.p + 1)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.x <= (model.p + 1)**2)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=(model.p + 1, model.x, model.p))
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=0 >= model.x - model.p)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.p >= model.x)
        self.assertTrue(model.c.upper is model.p)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.p >= model.x + 1)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=model.p + 1 >= model.x)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=(model.p + 1)**2 >= model.x)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=(None, model.x, model.p))
        self.assertTrue(model.c.upper is model.p)
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=(None, model.x + 1, model.p))
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=(None, model.x, model.p + 1))
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

        model.c = Constraint(expr=(1, model.x, model.p))
        self.assertEqual(model.c.equality, False)
        model.del_component(model.c)

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
    def test_nonlinear(self):
        m = ConcreteModel()
        m.x = Var()
        m.y = Var(initialize=0)

        m.c = Constraint(expr=m.x**2 == 16)
        m.x.set_value(1.0) # set an initial value
        calculate_variable_from_constraint(m.x, m.c, linesearch=False)
        self.assertAlmostEqual(value(m.x), 4)

        # test that infeasible constraint throws error
        m.d = Constraint(expr=m.x**2 == -1)
        m.x.set_value(1.25) # set the initial value
        with self.assertRaisesRegexp(
               RuntimeError, 'Iteration limit \(10\) reached'):
           calculate_variable_from_constraint(
               m.x, m.d, iterlim=10, linesearch=False)

        # same problem should throw a linesearch error if linesearch is on
        m.x.set_value(1.25) # set the initial value
        with self.assertRaisesRegexp(
                RuntimeError, "Linesearch iteration limit reached"):
           calculate_variable_from_constraint(
               m.x, m.d, iterlim=10, linesearch=True)

        # same problem should raise an error if initialized at 0
        m.x = 0
        with self.assertRaisesRegexp(
                RuntimeError, "Initial value for variable results in a "
                "derivative value that is very close to zero."):
            calculate_variable_from_constraint(m.x, m.c)

        # same problem should raise a value error if we are asked to
        # solve for a variable that is not present
        with self.assertRaisesRegexp(
                ValueError, "Variable derivative == 0"):
            calculate_variable_from_constraint(m.y, m.c)


        # should succeed with or without a linesearch
        m.e = Constraint(expr=(m.x - 2.0)**2 - 1 == 0)
        m.x.set_value(3.1)
        calculate_variable_from_constraint(m.x, m.e, linesearch=False)
        self.assertAlmostEqual(value(m.x), 3)

        m.x.set_value(3.1)
        calculate_variable_from_constraint(m.x, m.e, linesearch=True)
        self.assertAlmostEqual(value(m.x), 3)


        # we expect this to succeed with the linesearch
        m.f = Constraint(expr=1.0/(1.0+exp(-m.x))-0.5 == 0)
        m.x.set_value(3.0)
        calculate_variable_from_constraint(m.x, m.f, linesearch=True)
        self.assertAlmostEqual(value(m.x), 0)

        # we expect this to fail without a linesearch
        m.x.set_value(3.0)
        with self.assertRaisesRegexp(
                RuntimeError, "Newton's method encountered a derivative "
                "that was too close to zero"):
            calculate_variable_from_constraint(m.x, m.f, linesearch=False)

        # Calculate the bubble point of Benzene.  THe first step
        # computed by calculate_variable_from_constraint will make the
        # second term become complex, and the evaluation will fail.
        # This tests that the algorithm cleanly continues
        m = ConcreteModel()
        m.x = Var()
        m.pc = 48.9e5
        m.tc = 562.2
        m.psc = {'A': -6.98273,
                 'B': 1.33213,
                 'C': -2.62863,
                 'D': -3.33399,
        }
        m.p = 101325
        @m.Constraint()
        def f(m):
            return m.pc * \
                exp((m.psc['A'] * (1 - m.x / m.tc) +
                     m.psc['B'] * (1 - m.x / m.tc)**1.5 +
                     m.psc['C'] * (1 - m.x / m.tc)**3 +
                     m.psc['D'] * (1 - m.x / m.tc)**6
                 ) / (1 - (1 - m.x / m.tc))) - m.p == 0
        m.x.set_value(298.15)
        calculate_variable_from_constraint(m.x, m.f, linesearch=False)
        self.assertAlmostEqual(value(m.x), 353.31855602)
        m.x.set_value(298.15)
        calculate_variable_from_constraint(m.x, m.f, linesearch=True)
        self.assertAlmostEqual(value(m.x), 353.31855602)

        # Starting with an invalid guess (above TC) should raise an
        # exception
        m.x.set_value(600)
        output = six.StringIO()
        with LoggingIntercept(output, 'pyomo', logging.WARNING):
            if six.PY2:
                expectedException = ValueError
            else:
#.........这里部分代码省略.........

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
    def __init__(self, results,  threshold=None, k=10, solver="glpk", verbosity=0):
        """
        :param result: Epitope prediction result object from which the epitope selection should be performed
        :type result: :class:`~Fred2.Core.Result.EpitopePredictionResult`
        :param dict(str,float) threshold: A dictionary scoring the binding thresholds for each HLA
                                          :class:`~Fred2.Core.Allele.Allele` key = allele name; value = the threshold
        :param int k: The number of epitopes to select
        :param str solver: The solver to be used (default glpk)
        :param int verbosity: Integer defining whether additional debugg prints are made >0 => debug mode
        """

        #check input data
        if not isinstance(results, EpitopePredictionResult):
            raise ValueError("first input parameter is not of type EpitopePredictionResult")

        _alleles = copy.deepcopy(results.columns.values.tolist())

        #test if allele prob is set, if not set allele prob uniform
        #if only partly set infer missing values (assuming uniformity of missing value)
        prob = []
        no_prob = []
        for a in _alleles:
            if a.prob is None:
                no_prob.append(a)
            else:
                prob.append(a)

        if len(no_prob) > 0:
            #group by locus
            no_prob_grouped = {}
            prob_grouped = {}
            for a in no_prob:
                no_prob_grouped.setdefault(a.locus, []).append(a)
            for a in prob:
                prob_grouped.setdefault(a.locus, []).append(a)

            for g, v in no_prob_grouped.iteritems():
                total_loc_a = len(v)
                if g in prob_grouped:
                    remaining_mass = 1.0 - sum(a.prob for a in prob_grouped[g])
                    for a in v:
                        a.prob = remaining_mass/total_loc_a
                else:
                    for a in v:
                        a.prob = 1.0/total_loc_a
        probs = {a.name:a.prob for a in _alleles}
        if verbosity:
            for a in _alleles:
                print a.name, a.prob

        #start constructing model
        self.__solver = SolverFactory(solver)
        self.__verbosity = verbosity
        self.__changed = True
        self.__alleleProb = _alleles
        self.__k = k
        self.__result = None
        self.__thresh = {} if threshold is None else threshold

        # Variable, Set and Parameter preparation
        alleles_I = {}
        variations = []
        epi_var = {}
        imm = {}
        peps = {}
        cons = {}

        #unstack multiindex df to get normal df based on first prediction method
        #and filter for binding epitopes
        method = results.index.values[0][1]
        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
#.........这里部分代码省略.........

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
    def test_get_check_units_on_all_expressions(self):
        # this method is going to test all the expression types that should work
        # to be defensive, we will also test that we actually have the expected expression type
        # therefore, if the expression system changes and we get a different expression type,
        # we will know we need to change these tests

        uc = units
        kg = uc.kg
        m = uc.m

        model = ConcreteModel()
        model.x = Var()
        model.y = Var()
        model.z = Var()
        model.p = Param(initialize=42.0, mutable=True)

        # test equality
        self._get_check_units_ok(3.0*kg == 1.0*kg, uc, 'kg', expr.EqualityExpression)
        self._get_check_units_fail(3.0*kg == 2.0*m, uc, expr.EqualityExpression)

        # test inequality
        self._get_check_units_ok(3.0*kg <= 1.0*kg, uc, 'kg', expr.InequalityExpression)
        self._get_check_units_fail(3.0*kg <= 2.0*m, uc, expr.InequalityExpression)
        self._get_check_units_ok(3.0*kg >= 1.0*kg, uc, 'kg', expr.InequalityExpression)
        self._get_check_units_fail(3.0*kg >= 2.0*m, uc, expr.InequalityExpression)

        # test RangedExpression
        self._get_check_units_ok(inequality(3.0*kg, 4.0*kg, 5.0*kg), uc, 'kg', expr.RangedExpression)
        self._get_check_units_fail(inequality(3.0*m, 4.0*kg, 5.0*kg), uc, expr.RangedExpression)
        self._get_check_units_fail(inequality(3.0*kg, 4.0*m, 5.0*kg), uc, expr.RangedExpression)
        self._get_check_units_fail(inequality(3.0*kg, 4.0*kg, 5.0*m), uc, expr.RangedExpression)

        # test SumExpression, NPV_SumExpression
        self._get_check_units_ok(3.0*model.x*kg + 1.0*model.y*kg + 3.65*model.z*kg, uc, 'kg', expr.SumExpression)
        self._get_check_units_fail(3.0*model.x*kg + 1.0*model.y*m + 3.65*model.z*kg, uc, expr.SumExpression)

        self._get_check_units_ok(3.0*kg + 1.0*kg + 2.0*kg, uc, 'kg', expr.NPV_SumExpression)
        self._get_check_units_fail(3.0*kg + 1.0*kg + 2.0*m, uc, expr.NPV_SumExpression)

        # test ProductExpression, NPV_ProductExpression
        self._get_check_units_ok(model.x*kg * model.y*m, uc, 'kg * m', expr.ProductExpression)
        self._get_check_units_ok(3.0*kg * 1.0*m, uc, 'kg * m', expr.NPV_ProductExpression)
        self._get_check_units_ok(3.0*kg*m, uc, 'kg * m', expr.NPV_ProductExpression)
        # I don't think that there are combinations that can "fail" for products

        # test MonomialTermExpression
        self._get_check_units_ok(model.x*kg, uc, 'kg', expr.MonomialTermExpression)

        # test ReciprocalExpression, NPV_ReciprocalExpression
        self._get_check_units_ok(1.0/(model.x*kg), uc, '1 / kg', expr.ReciprocalExpression)
        self._get_check_units_ok(1.0/kg, uc, '1 / kg', expr.NPV_ReciprocalExpression)
        # I don't think that there are combinations that can "fail" for products

        # test PowExpression, NPV_PowExpression
        # ToDo: fix the str representation to combine the powers or the expression system
        self._get_check_units_ok((model.x*kg**2)**3, uc, 'kg ** 6', expr.PowExpression) # would want this to be kg**6
        self._get_check_units_fail(kg**model.x, uc, expr.PowExpression, UnitsError)
        self._get_check_units_fail(model.x**kg, uc, expr.PowExpression, UnitsError)
        self._get_check_units_ok(kg**2, uc, 'kg ** 2', expr.NPV_PowExpression)
        self._get_check_units_fail(3.0**kg, uc, expr.NPV_PowExpression, UnitsError)

        # test NegationExpression, NPV_NegationExpression
        self._get_check_units_ok(-(kg*model.x*model.y), uc, 'kg', expr.NegationExpression)
        self._get_check_units_ok(-kg, uc, 'kg', expr.NPV_NegationExpression)
        # don't think there are combinations that fan "fail" for negation

        # test AbsExpression, NPV_AbsExpression
        self._get_check_units_ok(abs(kg*model.x), uc, 'kg', expr.AbsExpression)
        self._get_check_units_ok(abs(kg), uc, 'kg', expr.NPV_AbsExpression)
        # don't think there are combinations that fan "fail" for abs

        # test the different UnaryFunctionExpression / NPV_UnaryFunctionExpression types
        # log
        self._get_check_units_ok(log(3.0*model.x), uc, None, expr.UnaryFunctionExpression)
        self._get_check_units_fail(log(3.0*kg*model.x), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_ok(log(3.0*model.p), uc, None, expr.NPV_UnaryFunctionExpression)
        self._get_check_units_fail(log(3.0*kg), uc, expr.NPV_UnaryFunctionExpression, UnitsError)
        # log10
        self._get_check_units_ok(log10(3.0*model.x), uc, None, expr.UnaryFunctionExpression)
        self._get_check_units_fail(log10(3.0*kg*model.x), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_ok(log10(3.0*model.p), uc, None, expr.NPV_UnaryFunctionExpression)
        self._get_check_units_fail(log10(3.0*kg), uc, expr.NPV_UnaryFunctionExpression, UnitsError)
        # sin
        self._get_check_units_ok(sin(3.0*model.x*uc.radians), uc, None, expr.UnaryFunctionExpression)
        self._get_check_units_fail(sin(3.0*kg*model.x), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_fail(sin(3.0*kg*model.x*uc.kg), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_ok(sin(3.0*model.p*uc.radians), uc, None, expr.NPV_UnaryFunctionExpression)
        self._get_check_units_fail(sin(3.0*kg), uc, expr.NPV_UnaryFunctionExpression, UnitsError)
        # cos
        self._get_check_units_ok(cos(3.0*model.x*uc.radians), uc, None, expr.UnaryFunctionExpression)
        self._get_check_units_fail(cos(3.0*kg*model.x), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_fail(cos(3.0*kg*model.x*uc.kg), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_ok(cos(3.0*model.p*uc.radians), uc, None, expr.NPV_UnaryFunctionExpression)
        self._get_check_units_fail(cos(3.0*kg), uc, expr.NPV_UnaryFunctionExpression, UnitsError)
        # tan
        self._get_check_units_ok(tan(3.0*model.x*uc.radians), uc, None, expr.UnaryFunctionExpression)
        self._get_check_units_fail(tan(3.0*kg*model.x), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_fail(tan(3.0*kg*model.x*uc.kg), uc, expr.UnaryFunctionExpression, UnitsError)
        self._get_check_units_ok(tan(3.0*model.p*uc.radians), uc, None, expr.NPV_UnaryFunctionExpression)
        self._get_check_units_fail(tan(3.0*kg), uc, expr.NPV_UnaryFunctionExpression, UnitsError)
#.........这里部分代码省略.........

# 需要导入模块: from pyomo.environ import ConcreteModel [as 别名]
# 或者: from pyomo.environ.ConcreteModel import p [as 别名]
  def __init__(self, results, threshold=None, dist_threshold=1.0, distance={}, expression={}, uncertainty={}, overlap=0, k=10, k_taa=0,
               solver="glpk", verbosity=0, include=[]):
        """
        :param results: Epitope prediction result object from which the epitope selection should be performed
        :type results: :class:`~Fred2.Core.Result.EpitopePredictionResult`
        :param dict(str,float) threshold: A dictionary scoring the binding thresholds for each HLA
                                          :class:`~Fred2.Core.Allele.Allele` key = allele name; value = the threshold
        :param float dist_threshold: Distance threshold: an epitope gets excluded if an epitope has dist-2-self score
                                    smaller or equal to this threshold for any HLA allele
        :param dict((str,str),float) distance: A dictionary with key: (peptide sequence, HLA name)
                                               and value the distance2self
        :param dict(str, float) expression: A dictionary with key: gene ID, and value: Gene expression
                                            in FPKM/RPKM or TPM
        :param dict((str,str),float) uncertainty: A dictionary with key (peptide seq, HLA name), and value the
                                                  associated uncertainty of the immunogenicity prediction
        :param int k: The number of epitopes to select
        :param int k_taa: The number of TAA epitopes to select
        :param str solver: The solver to be used (default glpk)
        :param int verbosity: Integer defining whether additional debug prints are made >0 => debug mode
        """

        # check input data
        if not isinstance(results, EpitopePredictionResult):
            raise ValueError("first input parameter is not of type EpitopePredictionResult")

        _alleles = results.columns.values.tolist()

        # generate abundance dictionary of HLA alleles default is 2.0 as values will be log2 transformed
        probs = {a.name:2.0 if a.get_metadata("abundance", only_first=True) is None else
                 a.get_metadata("abundance", only_first=True) for a in _alleles}

        # start constructing model
        self.__solver = SolverFactory(solver)
        self.__verbosity = verbosity
        self.__changed = True
        self.__alleleProb = _alleles
        self.__k = k
        self.__k_taa = k_taa
        self.__result = None
        self.__thresh = {} if threshold is None else threshold
        self.__included = include
        self.overlap=overlap

        # variable, set and parameter preparation
        alleles_I = {}
        variations = []
        epi_var = {}
        imm = {}
        peps = {}
        taa = []
        var_epi = {}
        cons = {}

        for a in _alleles:
            alleles_I.setdefault(a.name, set())

        # unstack multiindex df to get normal df based on first prediction method
        # and filter for binding epitopes
        method = results.index.values[0][1]
        res_df = results.xs(results.index.values[0][1], level="Method")

        # if predcitions are not available for peptides/alleles, replace by 0
        res_df.fillna(0, inplace=True)

        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)]

        res_df.fillna(0, inplace=True)
        # transform scores to 1-log50k(IC50) scores if neccassary
        # and generate mapping dictionaries for Set definitions
        for tup in res_df.itertuples():
            p = tup[0]
            seq = str(p)

            if any(distance.get((seq, a.name), 1.0) <= dist_threshold for a in _alleles):
                continue
            peps[seq] = p
            if p.get_metadata("taa",only_first=True):
                taa.append(seq)
            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)
#.........这里部分代码省略.........