Run GapFill command

def run(self):
        """Run GapFill command"""

        # Load compound information
        def compound_name(id):
            if id not in self._model.compounds:
                return id
            return self._model.compounds[id].properties.get('name', id)

        # Calculate penalty if penalty file exists
        penalties = {}
        if self._args.penalty is not None:
            for line in self._args.penalty:
                line, _, comment = line.partition('#')
                line = line.strip()
                if line == '':
                    continue
                rxnid, penalty = line.split(None, 1)
                penalties[rxnid] = float(penalty)

        core = set(self._mm.reactions)

        solver = self._get_solver(integer=True)
        default_comp = self._model.default_compartment
        epsilon = self._args.epsilon
        v_max = float(self._model.default_flux_limit)

        blocked = set()
        for compound in self._args.compound:
            if compound.compartment is None:
                compound = compound.in_compartment(default_comp)
            blocked.add(compound)

        if len(blocked) > 0:
            logger.info('Unblocking compounds: {}...'.format(
                ', '.join(text_type(c) for c in sorted(blocked))))
        else:
            logger.info(
                'Unblocking all compounds in model. Use --compound option to'
                ' unblock specific compounds.')
            blocked = set(self._mm.compounds)

        exclude = set()
        if self._model.biomass_reaction is not None:
            exclude.add(self._model.biomass_reaction)

        # Add exchange and transport reactions to database
        model_complete, weights = create_extended_model(
            self._model,
            db_penalty=self._args.db_penalty,
            ex_penalty=self._args.ex_penalty,
            tp_penalty=self._args.tp_penalty,
            penalties=penalties)

        implicit_sinks = not self._args.no_implicit_sinks

        logger.info('Searching for reactions to fill gaps')
        try:
            added_reactions, no_bounds_reactions = gapfill(
                model_complete, core, blocked, exclude, solver=solver,
                epsilon=epsilon, v_max=v_max, weights=weights,
                implicit_sinks=implicit_sinks,
                allow_bounds_expansion=self._args.allow_bounds_expansion)
        except GapFillError as e:
            self._log_epsilon_and_fail(epsilon, e)

        for reaction_id in sorted(self._mm.reactions):
            rx = self._mm.get_reaction(reaction_id)
            rxt = rx.translated_compounds(compound_name)
            print('{}\t{}\t{}\t{}'.format(reaction_id, 'Model', 0, rxt))

        for rxnid in sorted(added_reactions):
            rx = model_complete.get_reaction(rxnid)
            rxt = rx.translated_compounds(compound_name)
            print('{}\t{}\t{}\t{}'.format(
                rxnid, 'Add', weights.get(rxnid, 1), rxt))

        for rxnid in sorted(no_bounds_reactions):
            rx = model_complete.get_reaction(rxnid)
            rxt = rx.translated_compounds(compound_name)
            print('{}\t{}\t{}\t{}'.format(
                rxnid, 'Remove bounds', weights.get(rxnid, 1), rxt))