Convert exchange reactions in model to exchange compounds.

    Only exchange reactions in the extracellular compartment are converted.
    The extracelluar compartment must be defined for the model.

    Args:
        model: :class:`NativeModel`.
331934
def convert_exchange_to_compounds(model):
    
    # Build set of exchange reactions
    exchanges = set()
    for reaction in model.reactions:
        equation = reaction.properties.get('equation')
        if equation is None:
            continue

        if len(equation.compounds) != 1:
            # Provide warning for exchange reactions with more than
            # one compound, they won't be put into the exchange definition
            if (len(equation.left) == 0) != (len(equation.right) == 0):
                logger.warning('Exchange reaction {} has more than one'
                               ' compound, it was not converted to'
                               ' exchange compound'.format(reaction.id))
            continue

        exchanges.add(reaction.id)

    # Convert exchange reactions into exchange compounds
    for reaction_id in exchanges:
        equation = model.reactions[reaction_id].equation
        compound, value = equation.compounds[0]
        if compound.compartment != model.extracellular_compartment:
            continue

        if compound in model.exchange:
            logger.warning(
                'Compound {} is already defined in the exchange'
                ' definition'.format(compound))
            continue

        # We multiply the flux bounds by value in order to create equivalent
        # exchange reactions with stoichiometric value of one. If the flux
        # bounds are not set but the reaction is unidirectional, the implicit
        # flux bounds must be used.
        lower_flux, upper_flux = None, None
        if reaction_id in model.limits:
            _, lower, upper = model.limits[reaction_id]
            if lower is not None:
                lower_flux = lower * abs(value)
            if upper is not None:
                upper_flux = upper * abs(value)

        if lower_flux is None and equation.direction == Direction.Forward:
            lower_flux = 0
        if upper_flux is None and equation.direction == Direction.Reverse:
            upper_flux = 0

        # If the stoichiometric value of the reaction is reversed, the flux
        # limits must be flipped.
        if value > 0:
            lower_flux, upper_flux = (
                -upper_flux if upper_flux is not None else None,
                -lower_flux if lower_flux is not None else None)

        model.exchange[compound] = (
            compound, reaction_id, lower_flux, upper_flux)

        model.reactions.discard(reaction_id)
        model.limits.pop(reaction_id, None)
657981
def find_oxygen_reactions(model):
    """"""
    o2_in_model = helpers.find_met_in_model(model, "MNXM4")
    return set([rxn for met in model.metabolites for
                rxn in met.reactions if met.formula == "O2" or
                met in o2_in_model])
194616
def reciprocal_lattice_from_outcar( filename ): # from https://github.com/MaterialsDiscovery/PyChemia
    
    outcar = open(filename, "r").read()
    # just keeping the last component
    recLat = re.findall(r"reciprocal\s*lattice\s*vectors\s*([-.\s\d]*)",
                        outcar)[-1]
    recLat = recLat.split()
    recLat = np.array(recLat, dtype=float)
    # up to now I have, both direct and rec. lattices (3+3=6 columns)
    recLat.shape = (3, 6)
    recLat = recLat[:, 3:]
    return recLat
1050589
def parse_reaction_equation(equation_def, default_compartment):
    

    def parse_compound_list(l, compartment):
        """Parse a list of reactants or metabolites"""
        for compound_def in l:
            compound_id = compound_def.get('id')
            _check_id(compound_id, 'Compound')

            value = compound_def.get('value')
            if value is None:
                raise ParseError('Missing value for compound {}'.format(
                    compound_id))

            compound_compartment = compound_def.get('compartment')
            if compound_compartment is None:
                compound_compartment = compartment

            compound = Compound(compound_id, compartment=compound_compartment)
            yield compound, value

    if isinstance(equation_def, string_types):
        return parse_reaction_equation_string(
            equation_def, default_compartment)
    else:
        compartment = equation_def.get('compartment', default_compartment)
        reversible = bool(equation_def.get('reversible', True))
        left = equation_def.get('left', [])
        right = equation_def.get('right', [])
        if len(left) == 0 and len(right) == 0:
            raise ParseError('Reaction values are missing')

        return Reaction(Direction.Both if reversible else Direction.Forward,
                        parse_compound_list(left, compartment),
                        parse_compound_list(right, compartment))
900732
def translate_sbml_reaction(entry, new_id, compartment_map, compound_map):
    """"""
    new_entry = DictReactionEntry(entry, id=new_id)

    # Convert compound IDs in reaction equation
    if new_entry.equation is not None:
        compounds = []
        for compound, value in new_entry.equation.compounds:
            # Translate compartment to new ID, if available.
            compartment = compartment_map.get(
                compound.compartment, compound.compartment)
            new_compound = compound.translate(
                lambda name: compound_map.get(name, name)).in_compartment(
                    compartment)
            compounds.append((new_compound, value))

        new_entry.equation = Reaction(
            new_entry.equation.direction, compounds)

    # Get XHTML notes properties
    for key, value in iteritems(parse_xhtml_reaction_notes(entry)):
        if key not in new_entry.properties:
            new_entry.properties[key] = value

    return new_entry
54275
public static double calculateAverageBondLength(IChemModel model) {

        // empty models have to have a scale
        IAtomContainerSet moleculeSet = model.getMoleculeSet();
        if (moleculeSet == null) {
            IReactionSet reactionSet = model.getReactionSet();
            if (reactionSet != null) {
                return calculateAverageBondLength(reactionSet);
            }
            return 0.0;
        }

        return calculateAverageBondLength(moleculeSet);
    }
1363201
def internal_model_from_file(
            self, file_name, encoding='utf-8', debug=None,
            pre_ref_resolution_callback=None, is_main_model=True):
        
        model = None
        callback = pre_ref_resolution_callback

        if hasattr(self, "_tx_model_repository"):
            # metamodel has a global repo
            if not callback:
                def _pre_ref_resolution_callback(other_model):
                    from textx.scoping import GlobalModelRepository
                    filename = other_model._tx_filename
                    assert filename
                    # print("METAMODEL PRE-CALLBACK{}".format(filename))
                    other_model._tx_model_repository = GlobalModelRepository(
                        self._tx_model_repository.all_models)
                    self._tx_model_repository.all_models\
                        .filename_to_model[filename] = other_model

                callback = _pre_ref_resolution_callback
            if self._tx_model_repository.all_models.has_model(file_name):
                model = self._tx_model_repository.all_models\
                    .filename_to_model[file_name]

        if not model:
            # model not present (from global repo) -> load it
            model = self._parser_blueprint.clone().get_model_from_file(
                file_name, encoding, debug=debug,
                pre_ref_resolution_callback=callback,
                is_main_model=is_main_model)

        for p in self._model_processors:
            p(model, self)
        return model
1217186
def EXCHANGE(classical_reg1, classical_reg2):
    
    left = unpack_classical_reg(classical_reg1)
    right = unpack_classical_reg(classical_reg2)
    return ClassicalExchange(left, right)
1328630
def _assemble_conversion(self, stmt):
        """"""
        pybel_lists = ([], [])
        for pybel_list, agent_list in \
                            zip(pybel_lists, (stmt.obj_from, stmt.obj_to)):
            for agent in agent_list:
                node = _get_agent_grounding(agent)
                # TODO check for missing grounding?
                pybel_list.append(node)

        rxn_node_data = reaction(
            reactants=pybel_lists[0],
            products=pybel_lists[1],
        )
        obj_node = self.model.add_node_from_data(rxn_node_data)
        obj_edge = None  # TODO: Any edge information possible here?
        # Add node for controller, if there is one
        if stmt.subj is not None:
            subj_attr, subj_edge = _get_agent_node(stmt.subj)
            subj_node = self.model.add_node_from_data(subj_attr)
            edge_data_list = _combine_edge_data(pc.DIRECTLY_INCREASES,
                                           subj_edge, obj_edge, stmt.evidence)
            for edge_data in edge_data_list:
                self.model.add_edge(subj_node, obj_node, **edge_data)
1210858
def model_exchange(model):
    """"""
    # Determine the default flux limits. If the value is already at the
    # default it does not need to be included in the output.
    lower_default, upper_default = None, None
    if model.default_flux_limit is not None:
        lower_default = -model.default_flux_limit
        upper_default = model.default_flux_limit

    compounds = []
    for compound, reaction_id, lower, upper in sorted(
            itervalues(model.exchange)):
        d = OrderedDict([('id', compound.name)])
        if reaction_id is not None:
            d['reaction'] = reaction_id

        lower = _get_output_limit(lower, lower_default)
        upper = _get_output_limit(upper, upper_default)
        d.update(_generate_limit_items(lower, upper))

        compounds.append(d)

    return OrderedDict([('compounds', compounds)])
427406
def infer_compartment_entries(model):
    """"""
    compartment_ids = set()
    for reaction in model.reactions:
        equation = reaction.equation
        if equation is None:
            continue

        for compound, _ in equation.compounds:
            compartment = compound.compartment
            if compartment is None:
                compartment = model.default_compartment

            if compartment is not None:
                compartment_ids.add(compartment)

    for compartment in compartment_ids:
        if compartment in model.compartments:
            continue

        entry = DictCompartmentEntry(dict(id=compartment))
        model.compartments.add_entry(entry)
11925
def model_reaction_limits(model):
    """"""
    for reaction in sorted(model.reactions, key=lambda r: r.id):
        equation = reaction.properties.get('equation')
        if equation is None:
            continue

        # Determine the default flux limits. If the value is already at the
        # default it does not need to be included in the output.
        lower_default, upper_default = None, None
        if model.default_flux_limit is not None:
            if equation.direction.reverse:
                lower_default = -model.default_flux_limit
            else:
                lower_default = 0.0

            if equation.direction.forward:
                upper_default = model.default_flux_limit
            else:
                upper_default = 0.0

        lower_flux, upper_flux = None, None
        if reaction.id in model.limits:
            _, lower, upper = model.limits[reaction.id]
            lower_flux = _get_output_limit(lower, lower_default)
            upper_flux = _get_output_limit(upper, upper_default)

        if lower_flux is not None or upper_flux is not None:
            d = OrderedDict([('reaction', reaction.id)])
            d.update(_generate_limit_items(lower_flux, upper_flux))

            yield d
594705
def convert(cycle, pyroleLike, usedPyroles):
    """"""

    bonds = cycle.bonds
    atoms = cycle.atoms
    initialBondStates = []
    initialAtomStates = []
    _usedPyroles = {}
    for bond in bonds:
        # store the initial states but assume the
        # bond is aromatic
        initialBondStates.append((bond, bond.symbol,
                                  bond.bondorder, bond.bondtype,
                                  bond.aromatic, bond.stereo))
        # XXX FIX ME
        # until we get proper conjugation, aromatic bond orders
        # are 1.5
        bond.reset(':', bond.bondorder, 4, bond.fixed, bond.stereo)
        
    aromatized = 1
    for atom in atoms:
        initialAtomStates.append((atom, atom.aromatic))
        atom.aromatic = 1

        nonhydrogens = atom.sumBondOrders() + atom.charge

        # look for the lowest valence where we don't
        #  have to change the charge of the atom to
        #  fill the valences
        for valence in atom.valences:
            neededHydrogens = int(valence - nonhydrogens)
            if neededHydrogens >= 0:
                break
        else:
            # we can't change the aromaticity and have correct
            # valence.
            #
            # there is one special case of a five membered
            #  ring and a pyrole nitrogen like atom we need
            #  to look for.
            if len(cycle) == 5 and pyroleLike.has_key(atom.handle):
                _usedPyroles[atom.handle] = 1
            else:
                # nope, the valences don't work out so
                # we can't aromatize
                aromatized = 0
                break

    # sanity check, this should be true because of the
    # canBeAromatic routine above
    assert len(_usedPyroles) <=1, "Too many used pyroles!"
    
    cycle.aromatic = aromatized
    if not aromatized:
        for bond, symbol, order, bondtype, aromatic, stereo in initialBondStates:
            bond.reset(symbol, order, bondtype, bond.fixed, stereo)

        for atom, aromatic in initialAtomStates:
            atom.aromatic = aromatic
    else:
        # we used some pyroles, we'll have to send these to
        # the valence checker later
        usedPyroles.update(_usedPyroles)

    return aromatized
457755
def load_model(cls, database, reaction_iter=None, exchange=None,
                   limits=None, v_max=None):
        

        model_args = {}
        if v_max is not None:
            model_args['v_max'] = v_max

        model = cls(database, **model_args)
        if reaction_iter is None:
            reaction_iter = iter(database.reactions)
        for reaction_id in reaction_iter:
            model.add_reaction(reaction_id)

        # Apply reaction limits
        if limits is not None:
            for reaction_id, lower, upper in limits:
                if lower is not None:
                    model.limits[reaction_id].lower = lower
                if upper is not None:
                    model.limits[reaction_id].upper = upper

        # TODO: Currently we just create a new exchange reaction in the
        # database and add it to the model. Ideally, we should not modify
        # the database. The exchange reaction could be created on the
        # fly when required.
        if exchange is not None:
            for compound, reaction_id, lower, upper in exchange:
                # Create exchange reaction
                if reaction_id is None:
                    reaction_id = create_exchange_id(
                        model.database.reactions, compound)
                model.database.set_reaction(
                    reaction_id, Reaction(Direction.Both, {compound: -1}))
                model.add_reaction(reaction_id)
                if lower is not None:
                    model.limits[reaction_id].lower = lower
                if upper is not None:
                    model.limits[reaction_id].upper = upper

        return model
254873
def apply_transformation(self, structure, return_ranked_list=False):
        
        comp = structure.composition
        logger.info("Composition: %s" % comp)

        for sp in comp:
            try:
                sp.oxi_state
            except AttributeError:
                analyzer = BVAnalyzer()
                structure = analyzer.get_oxi_state_decorated_structure(
                    structure)
                comp = structure.composition
                break

        ox = self.dopant.oxi_state
        radius = self.dopant.ionic_radius

        compatible_species = [
            sp for sp in comp if sp.oxi_state == ox and
            abs(sp.ionic_radius / radius - 1) < self.ionic_radius_tol]

        if (not compatible_species) and self.alio_tol:
            # We only consider aliovalent doping if there are no compatible
            # isovalent species.
            compatible_species = [
                sp for sp in comp
                if abs(sp.oxi_state - ox) <= self.alio_tol and
                abs(sp.ionic_radius / radius - 1) < self.ionic_radius_tol and
                sp.oxi_state * ox >= 0]

        if self.allowed_doping_species is not None:
            # Only keep allowed doping species.
            compatible_species = [
                sp for sp in compatible_species
                if sp in [get_el_sp(s) for s in self.allowed_doping_species]]

        logger.info("Compatible species: %s" % compatible_species)

        lengths = structure.lattice.abc
        scaling = [max(1, int(round(math.ceil(self.min_length / x))))
                   for x in lengths]
        logger.info("Lengths are %s" % str(lengths))
        logger.info("Scaling = %s" % str(scaling))

        all_structures = []
        t = EnumerateStructureTransformation(**self.kwargs)

        for sp in compatible_species:
            supercell = structure * scaling
            nsp = supercell.composition[sp]
            if sp.oxi_state == ox:
                supercell.replace_species({sp: {sp: (nsp - 1) / nsp,
                                                self.dopant: 1 / nsp}})
                logger.info("Doping %s for %s at level %.3f" % (
                    sp, self.dopant, 1 / nsp))
            elif self.codopant:
                codopant = _find_codopant(sp, 2 * sp.oxi_state - ox)
                supercell.replace_species({sp: {sp: (nsp - 2) / nsp,
                                                self.dopant: 1 / nsp,
                                                codopant: 1 / nsp}})
                logger.info("Doping %s for %s + %s at level %.3f" % (
                    sp, self.dopant, codopant, 1 / nsp))
            elif abs(sp.oxi_state) < abs(ox):
                # Strategy: replace the target species with a
                # combination of dopant and vacancy.
                # We will choose the lowest oxidation state species as a
                # vacancy compensation species as it is likely to be lower in
                # energy
                sp_to_remove = min([s for s in comp if s.oxi_state * ox > 0],
                                   key=lambda ss: abs(ss.oxi_state))

                if sp_to_remove == sp:
                    common_charge = lcm(int(abs(sp.oxi_state)), int(abs(ox)))
                    ndopant = common_charge / abs(ox)
                    nsp_to_remove = common_charge / abs(sp.oxi_state)
                    logger.info("Doping %d %s with %d %s." %
                                (nsp_to_remove, sp, ndopant, self.dopant))
                    supercell.replace_species(
                        {sp: {sp: (nsp - nsp_to_remove) / nsp,
                              self.dopant: ndopant / nsp}})
                else:
                    ox_diff = int(abs(round(sp.oxi_state - ox)))
                    vac_ox = int(abs(sp_to_remove.oxi_state))
                    common_charge = lcm(vac_ox, ox_diff)
                    ndopant = common_charge / ox_diff
                    nx_to_remove = common_charge / vac_ox
                    nx = supercell.composition[sp_to_remove]
                    logger.info("Doping %d %s with %s and removing %d %s." %
                                (ndopant, sp, self.dopant,
                                 nx_to_remove, sp_to_remove))
                    supercell.replace_species(
                        {sp: {sp: (nsp - ndopant) / nsp,
                              self.dopant: ndopant / nsp},
                         sp_to_remove: {
                             sp_to_remove: (nx - nx_to_remove) / nx}})
            elif abs(sp.oxi_state) > abs(ox):
                # Strategy: replace the target species with dopant and also
                # remove some opposite charged species for charge neutrality
                if ox > 0:
                    sp_to_remove = max(supercell.composition.keys(),
                                       key=lambda el: el.X)
                else:
                    sp_to_remove = min(supercell.composition.keys(),
                                       key=lambda el: el.X)
                # Confirm species are of opposite oxidation states.
                assert sp_to_remove.oxi_state * sp.oxi_state < 0

                ox_diff = int(abs(round(sp.oxi_state - ox)))
                anion_ox = int(abs(sp_to_remove.oxi_state))
                nx = supercell.composition[sp_to_remove]
                common_charge = lcm(anion_ox, ox_diff)
                ndopant = common_charge / ox_diff
                nx_to_remove = common_charge / anion_ox
                logger.info("Doping %d %s with %s and removing %d %s." %
                            (ndopant, sp, self.dopant,
                             nx_to_remove, sp_to_remove))
                supercell.replace_species(
                    {sp: {sp: (nsp - ndopant) / nsp,
                          self.dopant: ndopant / nsp},
                     sp_to_remove: {sp_to_remove: (nx - nx_to_remove) / nx}})

            ss = t.apply_transformation(
                supercell, return_ranked_list=self.max_structures_per_enum)
            logger.info("%s distinct structures" % len(ss))
            all_structures.extend(ss)

        logger.info("Total %s doped structures" % len(all_structures))
        if return_ranked_list:
            return all_structures[:return_ranked_list]

        return all_structures[0]["structure"]
107116
def _converter(self, value):
        
        if type(value) is not self.related_model_cls.Collection:
            value = self.related_model_cls.Collection([
                self._get_model_instance(self.related_model_cls, item)
                for item in value])
        return value
581679
def dissolved(self, concs):
        
        new_concs = concs.copy()
        for r in self.rxns:
            if r.has_precipitates(self.substances):
                net_stoich = np.asarray(r.net_stoich(self.substances))
                s_net, s_stoich, s_idx = r.precipitate_stoich(self.substances)
                new_concs -= new_concs[s_idx]/s_stoich * net_stoich
        return new_concs