""" methods for comparison of structures using rmsd """
# Contributors: Leif D. Jacobson
import abc
import collections
import hashlib
import itertools
import networkx as nx
import numpy as np
from networkx.algorithms import isomorphism
import schrodinger.application.jaguar.autots_constants as autots_constants
from schrodinger.application.jaguar import file_logger
from schrodinger.application.jaguar.autots_bonding import active_atom_pairs
from schrodinger.application.jaguar.autots_bonding import clean_st
from schrodinger.application.jaguar.autots_bonding import copy_autots_properties
from schrodinger.application.jaguar.autots_bonding import \
copy_autots_st_properties
from schrodinger.application.jaguar.autots_bonding import get_mmlewis_bonding
from schrodinger.application.jaguar.autots_bonding import remove_formal_charges
from schrodinger.application.jaguar.autots_exceptions import IsomerError
from schrodinger.application.jaguar.autots_exceptions import TemplateAtomMapFailure
from schrodinger.comparison import chirality
from schrodinger.comparison import are_conformers
from schrodinger.comparison import are_isomers
from schrodinger.comparison.atom_mapper import BaseAtomMapper
from schrodinger.comparison.atom_mapper import _register_file
from schrodinger.comparison.atom_mapper import _fill_map
from schrodinger.comparison.atom_mapper import ACTIVE_ATOM
from schrodinger.comparison.exceptions import ConformerError
from schrodinger.structutils import build
from schrodinger.structutils import rmsd
from schrodinger.structutils import transform
hashfunc = hashlib.sha256
MapScore = collections.namedtuple("MapScore", ["score", "atom_map"])
[docs]def align_all_molecules(st1, st2, sort_rms=True):
"""
Inspects whether or not all the molecules in reactant
are conformers of those in product.
If they are, each molecule in product is superimposed onto the
corresponding reactant molecule.
:type st1: schrodinger.structure.Structure instance
:param st1: structure 1
:type st2: schrodinger.structure.Structure instance
:param st2: structure 2
:type sort_rms: boolean
:param sort_rms: If True the structure in which the inter-molecule distances are
larger is superimposes onto the smaller, see sort_by_centroid_distance,
which does the sorting.
"""
if are_conformers(st1, st2, use_lewis_structure=False):
if sort_rms:
destination, source = sort_by_centroid_distance([st1, st2])
else:
source = st2
destination = st1
for mol in st1.molecule:
atlist = mol.getAtomIndices()
rmsd.superimpose(
destination, atlist, source, atlist, move_which=rmsd.MOLECULES)
[docs]def sort_by_centroid_distance(strs):
"""
Given a list of structures sort the list in non-decreasing order
using the rms distance between molecules in the structures.
e.g. for two structures which are each water dimers, the water dimer with the smaller
intermolecular distance would appear first in the list
:type strs: list of schrodinger.structure.Structure instances
:param strs: list of structures to be sorted
:rtype: list of schrodinger.structure.Structure instances
:return: sorted list
"""
nstrs = len(strs)
dists = [0.0 for i in range(nstrs)]
for ii, st in enumerate(strs):
nmol = len(st.molecule)
if nmol > 1:
mol_list = []
for mol in st.molecule:
mol_list.append(mol.extractStructure())
for i in range(nmol):
moli = mol_list[i]
xi, yi, zi, junk = transform.get_centroid(moli)
for j in range(i + 1, nmol):
molj = mol_list[j]
xj, yj, zj, junk = transform.get_centroid(molj)
dists[
ii] += (xi - xj)**2.0 + (yi - yj)**2.0 + (zi - zj)**2.0
dists[ii] /= float(nmol)
ind = np.argsort(dists)
out_strs = []
for i in ind:
out_strs.append(strs[i])
return out_strs
[docs]def align_path_strs(path_strs):
"""
Align the structures in a path.
:type path_strs: list of structures
:param path_strs: the structures along the path
"""
if not path_strs:
return
atlist = [i + 1 for i in range(path_strs[0].atom_total)]
# this imposes structure 1 on 0, 2 on 1, 3 on 2, etc.
for st1, st2 in zip(path_strs[1:], path_strs[:-1]):
rms = rmsd.superimpose(st2, atlist, st1, atlist, move_which=rmsd.CT)
[docs]def order_atoms(reactant, product, debug=False):
"""
Renumber atoms in reactant and product in a consistent fashion.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant structure
:type product: schrodinger.structure.Structure
:param product: reactant structure
:rtype: tuple of two structures
:return: renumbered reactant and product
"""
r_copy = reactant.copy()
p_copy = product.copy()
# first try using chirality
atlist = range(1, r_copy.atom_total + 1)
mapper = AutoTSAtomMapper(
optimize_mapping=True, use_chirality=True, debug=debug)
try:
r_new, p_new = mapper.reorder_structures(r_copy, atlist, p_copy, atlist)
except ConformerError as e:
mapper = AutoTSAtomMapper(
optimize_mapping=True, use_chirality=False, debug=debug)
r_new, p_new = mapper.reorder_structures(r_copy, atlist, p_copy, atlist)
copy_autots_st_properties(reactant, r_new)
copy_autots_st_properties(product, p_new)
return r_new, p_new
[docs]class AutoTSAtomMapper(BaseAtomMapper):
"""
Atom mapper used by AutoTS for reordering atoms of conformers.
Uses element types and optionally chirality to equate atoms
and all bonds are considered equivalent.
Stereochemistry on atoms which are marked with the property
%s are ignored when comparing stereochemistry.
This property is used to mark the atoms in the reaction center.
Usage example:
mapper = AutoTSAtomMapper(optimize_mapping=True, use_chirality=False)
st1, st2 = mapper.reorder_atoms(st1, range(1, st1.atom_total + 1), st2, range(st2.atom_total + 1))
""" % autots_constants.IGNORE_ATOM_STEREO
RMSD_THRESH = 1.0e-8
[docs] def __init__(self, optimize_mapping, use_chirality, debug=False):
"""
:type optimize_mapping: boolean
:param optimize_mapping: if True search over all bijections to
find the one with lowest score
:type use_chirality: boolean
:param use_chirality: if True, in addition to element type use
chirality (where defined) to equate atoms.
"""
super(AutoTSAtomMapper, self).__init__(optimize_mapping, debug)
self.use_chirality = use_chirality
[docs] def initialize_atom_types(self, st, invert_stereocenters=False):
"""
Initialize the atom types
:type st: Structure
:param st: structure containing atoms
:type invert_stereocenters: boolean
:param invert_stereocenters: whether or not R/S labels should be flipped
"""
def atom_label(atom):
"""
creates a label which is typically just the element name
but also contains a special atom class integer if the atom
has the property %s
""" % autots_constants.AUTOTS_ATOM_CLASS
atom_class = atom.property.get(
str(autots_constants.AUTOTS_ATOM_CLASS), "")
label = atom.element
if atom_class:
label += "-class=%s" % atom_class
return label
# uses element-chirality
if self.use_chirality:
ch = chirality.get_chirality(st)
if invert_stereocenters:
self.invert_chirality(ch)
for at in st.atom:
# have to get rid of these atom number chiralities here
# since not using Lewis structure, E/Z isomerism doesn't apply
if ch[at.index - 1] in ["ANR", "ANS", "E", "Z"]:
ch[at.index - 1] = None
self.set_atom_type(
at,
atom_label(at) + "-chiral=%s" % ch[at.index - 1])
else:
for at in st.atom:
self.set_atom_type(at, atom_label(at))
[docs] def score_mapping(self, st1, st2, atset):
"""
Score a mapping over the atoms listed in atlist.
This returns a number of chirality mis-alignments and an rms
mis-alignments are neglected for atoms which are marked with
atom property %s
:type st1: Structure
:param st1: first structure
:type st2: Structure
:param st2: second structure
:type atset: set of ints
:param atset: set of atoms that have been mapped (indexes refer to both structures)
:rtype: 3-tuple of (non-negative) int, (non-positive) int, and real
:return: the number of misaligned atom-numbering stereocenters,
the negative of the number of active atoms, and the rmsd
""" % autots_constants.IGNORE_ATOM_STEREO
atlist = list(atset)
# move a copy
st1_working = st1.copy()
st2_working = st2.copy()
rms = rmsd.superimpose_substructure_molecules(st1_working, atlist,
st2_working, atlist)
# according to cProfile this is the bottleneck
ch1 = chirality.get_atom_numbering_chirality(st1_working)
ch2 = chirality.get_atom_numbering_chirality(st2_working)
misalignments = 0
# only compare chirality for atoms which have all their
# neighbors mapped
atlist_neighbors_mapped = [
i for i in atlist
if all(
[at.index in atset for at in st1_working.atom[i].bonded_atoms])
]
for i in atlist_neighbors_mapped:
ignore_stereo = st1.atom[i].property.get(autots_constants.IGNORE_ATOM_STEREO, False) or \
st2.atom[i].property.get(autots_constants.IGNORE_ATOM_STEREO, False)
if (chirality.valid_chiral_comp(ch1[i - 1], ch2[i - 1]) and
ch1[i - 1] != ch2[i - 1]) and not ignore_stereo:
misalignments += 1
# In order to encourage larger n-membered TS's, we prefer solutions that
# increase the number of active atoms (even at the expense of RMSD).
# We include the negative of the number of active atoms to maximize this when
# we take the minimum score later.
# renumber.py has added this property to the active atoms for the purposes of
# renumbering. It also handles removing the label when renumbering is complete.
active_indices = set()
for i in atlist:
if st1.atom[i].property.get(autots_constants.IGNORE_ATOM_STEREO, False) or \
st2.atom[i].property.get(autots_constants.IGNORE_ATOM_STEREO, False):
active_indices.add(i)
score = (misalignments, -len(active_indices), rms)
if self.debug:
ignored_atoms = []
for i in atlist_neighbors_mapped:
if st1.atom[i].property.get(autots_constants.IGNORE_ATOM_STEREO, False) or \
st2.atom[i].property.get(autots_constants.IGNORE_ATOM_STEREO, False):
ignored_atoms.append(i)
print("ignoring stereo for the following atoms", ignored_atoms)
st1_working.append("scoring.mae")
st2_working.append("scoring.mae")
file_logger.register_file("scoring.mae")
print("score", score)
print("misalignments:")
for i in atlist_neighbors_mapped:
if ch1[i - 1] != ch2[i - 1]:
print(i, ch1[i - 1], ch2[i - 1])
return score
[docs] def score_is_equivalent(self, score1, score2):
"""
Here we declare 2 scores equivalent if the same number
of chirality mismatches are present and the rmsd
difference is within RMSD_THRESH. This resolves
machine dependent descrepancies in the chosen map.
:type score1: tuple
:param score1: the first score (chirality mismatches, active atoms, rms)
:type score2: tuple
:param score2: the first score (chirality mismatches, active atoms, rms)
:return boolean indicating if the score is the same or different
"""
same_chirality = score1[0] == score2[0]
same_active = score1[1] == score2[1]
same_rms = abs(score1[2] - score2[2]) < self.RMSD_THRESH
return same_chirality and same_rms and same_active
[docs]class AutoTSTemplateAtomMapper(AutoTSAtomMapper):
"""
Atom mapper used by AutoTS's template code for reordering atoms
to align template and input.
Uses element types to equate atoms and all bonds are considered equivalent.
Stereochemistry on atoms which are marked with the property
%s are ignored when comparing stereochemistry.
This differs from the AutoTSAtomMapper in that it simultaneously renumbers
4 structures instead of 2, assuming two pairs of those structures are already
consistently numbered. Moreover, the return value of this class is the optimal map
that accomplishes the renumbering, as opposed to the renumbered Structures themselves.
Usage example:
mapper = AutoTSTemplateAtomMapper(optimize_mapping=True, use_chirality=False)
map = mapper.choose_template_map(reactant, product, input_indexes, r_template, p_template, template_indexes)
""" % autots_constants.IGNORE_ATOM_STEREO
RMSD_THRESH = 1.0e-8
[docs] def initialize_atom_types(self, st, active_ats, invert_stereocenters=False):
"""
Initialize the atom types
:type st: Structure
:param st: structure containing atoms
:type active_ats: list of ints
:param active_ats: atom indices of the reaction center
:type invert_stereocenters: boolean
:param invert_stereocenters: whether or not R/S labels should be flipped
"""
for at in st.atom:
if autots_constants.AUTOTS_ATOM_CLASS in at.property:
del at.property[autots_constants.AUTOTS_ATOM_CLASS]
super().initialize_atom_types(st, invert_stereocenters)
for at in active_ats:
at_type = self.get_atom_type(st.atom[at])
at_type += ACTIVE_ATOM
self.set_atom_type(st.atom[at], at_type)
def _order_atoms(self, st1, atset1, st2, atset2, st3, st4, input_rc,
template_rc):
"""
Order atoms by first recursively removing
topologically degenerate terminal atoms
:type st1: Structure
:param st1: the first structure
:type atset1: set of ints
:param atset1: set of atom indexes defining the first substructure
:type st2: Structure
:param st2: the second structure
:type atset2: set of ints
:param atset2: set of atom indexes defining the second substructure
:type st3: Structure
:param st3: the third structure, already assumed to be numbered
consistently with st1
:type st4: Structure
:param st4: the fourth structure, already assumed to be numbered
consistently with st2
:type input_rc: set of ints
:param input_rc: reaction center atom indices of st1, st3
:type template_rc: set of ints
:param template_rc: reaction center atom indices of st2, st4
:return: the four structures with structure 2, 4 having had atoms reordered
and the overall map which accomplished this
"""
st1_renorm, k1, e1 = self._eliminate_degenerate_terminal_atoms(
st1, atset1)
st2_renorm, k2, e2 = self._eliminate_degenerate_terminal_atoms(
st2, atset2)
st3_renorm, k3, e3 = self._eliminate_degenerate_terminal_atoms(
st3, atset1)
st4_renorm, k4, e4 = self._eliminate_degenerate_terminal_atoms(
st4, atset2)
# check elimination for similarity or if we tried to eliminate the rxn center
if not (self.isomeric_atom_sets(st1, e1, st2, e2) and
self.isomeric_atom_sets(st2, e2, st3, e3) and
self.isomeric_atom_sets(st3, e3, st4, e4)) or len(
e1.intersection(input_rc)) > 0 or len(
e2.intersection(template_rc)) > 0:
e1.clear()
e2.clear()
e3.clear()
e4.clear()
# if we've eliminated any terminal atoms recurse
if e1 and e2 and e3 and e4:
st1_reord, st2_reord, st3_reord, st4_reord, overall_map = self._order_atoms(
st1_renorm, k1, st2_renorm, k2, st3_renorm, st4_renorm,
input_rc, template_rc)
# order eliminated atoms
# because st1-st3 (r and p) and st2-st4 (r_template and p_template) are consistently numbered
# one actually need only consider st1 and st2 and then apply the same map selection to the
# other two st.
st1_reord, st2_reord, st3_reord, st4_reord, app_map = self._append_terminal_atoms(
st1_reord, e1, st2_reord, e2, st3_reord, st4_reord, k1)
overall_map.update(app_map)
# order core
else:
st1_reord, st2_reord, st3_reord, st4_reord, overall_map = self._order_core_atoms(
st1, atset1, st2, atset2, st3, st4, input_rc, template_rc)
return st1_reord, st2_reord, st3_reord, st4_reord, overall_map
def _append_terminal_atoms(self, st1, term1, st2, term2, st3, st4, core):
"""
Map the terminal atoms in the list term1 and term2.
The atom indexes in the list core should be consistently
ordered between st1 and st2.
:type st1: Structure
:param st1: the first structure
:type term1: set of ints
:param term1: set of terminal atom indexes in st1
:type st2: Structure
:param st2: the second structure
:type term2: set of ints
:param term2: set of terminal atom indexes in st2
:type st3: Structure
:param st3: the third structure, already assumed to be numbered
consistently with st1
:type st4: Structure
:param st4: the fourth structure, already assumed to be numbered
consistently with st2
:type core: set of atoms in core structure which are already
numbered consistently
:return: four new structures, map
"""
# keeps track of which atoms have been mapped
overall_map = dict()
have_mapped = core.copy()
for icore in core:
types1 = self._group_neighbors(st1.atom[icore], term1)
types2 = self._group_neighbors(st2.atom[icore], term2)
# check that name of groups and lengths are the same
group_cnts1 = {grp: len(types1[grp]) for grp in types1}
group_cnts2 = {grp: len(types2[grp]) for grp in types2}
if group_cnts1 != group_cnts2:
raise ConformerError(
"Groups are not the same when appending terminal atoms")
# each group has a list of mappings
group_mappings = collections.defaultdict(list)
for group in types1:
for p in itertools.permutations(types2[group], len(
types2[group])):
group_mappings[group].append(list(zip(types1[group], p)))
# the set of atoms getting mapped
mapped_atoms = set()
for lst in types1.values():
mapped_atoms.update(lst)
# cartesian product of all groups to get
# all mappings of this atoms neighbors
best_score = None
best_maps = []
for mapping in itertools.product(*group_mappings.values()):
# get a mapping of terminal atoms attached to this core
full_list = []
for group_mapping in mapping:
full_list.extend(group_mapping)
# set initial map to None
full_map = [None for i in range(st2.atom_total)]
# map everything that has been mapped
for j in have_mapped:
full_map[j - 1] = j
# map this grouping of atoms
for group_mapping in mapping:
for i, j in group_mapping:
full_map[i - 1] = j
# fill in Nones
_fill_map(full_map)
# RMS is computed over core atom and added terminals
# whereas chirality need only be computed for the new core atom
# this is because chirality is determined by core + terminal
st2_reord = build.reorder_atoms(st2, full_map)
st4_reord = build.reorder_atoms(st4, full_map)
all_mapped = have_mapped.union(mapped_atoms)
score_r = self.score_mapping(st1, st2_reord, all_mapped)
score_p = self.score_mapping(st3, st4_reord, all_mapped)
score = tuple(map(sum, zip(score_r, score_p)))
map_score = MapScore(score, full_map)
self._update_map_list(best_maps, map_score)
if not self.optimize_mapping:
break
## now pick the best map. If there are ties, choose (entirely arbitrarily)
## the one that preserves the relative ordering of the atoms being reordered
## as best as possible
if best_maps:
best_score, best_mapping = min(
best_maps, key=lambda bij: bij.atom_map)
# save the overall map
for at in mapped_atoms:
overall_map[at] = st2.atom[best_mapping[at - 1]].property[
self.MAPPER_INDEX]
if self.debug:
if best_maps:
print("Best app maps")
for mapping in best_maps:
print(str(mapping.score))
print("overall map in append", overall_map)
st2 = build.reorder_atoms(st2, best_mapping)
st4 = build.reorder_atoms(st4, best_mapping)
# we just mapped terminal atoms
have_mapped.update(mapped_atoms)
return st1, st2, st3, st4, overall_map
def _order_core_atoms(self, st1, atset1, st2, atset2, st3, st4, input_rc,
template_rc):
"""
re-order the atoms atset2 in st2/st4 to be consistent
with atset1 in st1/st3. Graphs made from the atom
indexes listed in atset1 and atset2 must be isomorphic.
:type st1: Structure
:param st1: the first structure
:type atset1: set of ints
:param atset1: set of atom indexes defining the first substructure
:type st2: Structure
:param st2: the second structure
:type atset2: list of ints
:param atset2: list of atom indexes defining the second substructure
:type st3: Structure
:param st3: the third structure, already assumed to be numbered
consistently with st1
:type st4: Structure
:param st4: the fourth structure, already assumed to be numbered
consistently with st2
:type input_rc: set of ints
:param input_rc: reaction center atom indices of st1, st3
:type template_rc: set of ints
:param template_rc: reaction center atom indices of st2, st4
:return: the four structures with structure 2, 4 having had atoms reordered
and the overall map which accomplished this
"""
if len(atset1) != len(atset2):
raise ConformerError("length of atom lists in core region differ")
r_graph = self.st_to_graph(st1, atset1)
r_template_graph = self.st_to_graph(st2, atset2)
p_graph = self.st_to_graph(st3, atset1)
p_template_graph = self.st_to_graph(st4, atset2)
r_matcher = isomorphism.GraphMatcher(r_graph, r_template_graph,
self.comparator)
p_matcher = isomorphism.GraphMatcher(p_graph, p_template_graph,
self.comparator)
best_maps = []
overall_map = dict()
for trial_map in _template_map_iter(r_matcher, p_matcher, input_rc,
template_rc):
atom_map = [None for i in range(st2.atom_total)]
for k, v in trial_map.items():
atom_map[k - 1] = v
_fill_map(atom_map)
r_t_reord = build.reorder_atoms(st2, atom_map)
p_t_reord = build.reorder_atoms(st4, atom_map)
score_r = self.score_mapping(st1, r_t_reord, atset1)
score_p = self.score_mapping(st3, p_t_reord, atset1)
score = tuple(map(sum, zip(score_r, score_p)))
map_score = MapScore(score, atom_map)
self._update_map_list(best_maps, map_score)
if best_maps:
best_score, best_mapping = min(
best_maps, key=lambda bij: bij.atom_map)
else:
raise TemplateAtomMapFailure
r_t_reord = build.reorder_atoms(st2, best_mapping)
p_t_reord = build.reorder_atoms(st4, best_mapping)
for i in atset1:
overall_map[i] = best_mapping[i - 1]
if self.debug:
print("best_score", best_score)
print(overall_map)
return st1, r_t_reord, st3, p_t_reord, overall_map
[docs] def choose_template_map(self, reactant, product, input_indexes, r_template,
p_template, template_indexes):
"""
Determine the optimal map from the input reactant and product structures
to the template reactant and product structures.
:type reactant: Structure instance
:param reactant: input reactant structure
:type product: Structure instance
:param product: input product structure
:type input_indexes: RxnIndxDecomp instance
:param input_indexes: breakdown of input atom indexes into core,
reaction center
:type r_template: Structure instance
:param r_template: template reactant structure
:type p_template: Structure instance
:param p_template: template product structure
:type template_indexes: RxnIndxDecomp instance
:param template_indexes: breakdown of template atom indexes into core,
reaction center
:rtype: dict
:return: mapping from input number to template numbering
"""
r_working = reactant.copy()
p_working = product.copy()
r_t_working = r_template.copy()
p_t_working = p_template.copy()
for sts, indexes in zip(
[[r_working, p_working], [r_t_working, p_t_working]],
[input_indexes, template_indexes]):
for st in sts:
self.initialize_atom_types(st, indexes.rxn_center)
for at in indexes.rxn_center:
st.atom[at].property[
autots_constants.IGNORE_ATOM_STEREO] = True
for at in st.atom:
at.property[self.MAPPER_INDEX] = at.index
# Because our scoring works by reordering the structures, the actual map
# determination must involve the system with more atoms being mapped onto the
# system with fewer atoms (since you can't reorder the system with fewer atoms to
# have atom indices greater than its number of atoms.
# We simply reverse the map if need be to ensure a map from input to template
if r_working.atom_total <= r_t_working.atom_total:
r_working, r_t_working, p_working, p_t_working, overall_map = self._order_atoms(
r_working, set(input_indexes.core), r_t_working,
set(template_indexes.core), p_working, p_t_working,
set(input_indexes.rxn_center), set(template_indexes.rxn_center))
if self.debug:
r_t_working.write("rt_reord.mae")
file_logger.register_file("rt_reord.mae")
p_t_working.write("pt_reord.mae")
file_logger.register_file("pt_reord.mae")
else:
r_t_working, r_working, p_t_working, p_working, reverse_map = self._order_atoms(
r_t_working, set(template_indexes.core), r_working,
set(input_indexes.core), p_t_working, p_working,
set(template_indexes.rxn_center), set(input_indexes.rxn_center))
overall_map = {v: k for k, v in reverse_map.items()}
if self.debug:
r_working.write("rt_reord_rev.mae")
file_logger.register_file("rt_reord_rev.mae")
p_working.write("pt_reord_rev.mae")
file_logger.register_file("pt_reord_rev.mae")
if self.debug:
print("final overall_map", overall_map)
return overall_map
def _template_map_iter(r_matcher, p_matcher, input_rc, template_rc):
"""
Iterator over valid template isomorphism atom maps.
Valid atom maps must exist between both reactants and product and map
the reaction center to the reaction center.
Since template r and p are consistently numbered and input r and p are consistently
numbered, only maps that appear in both lists are good maps to use.
That is, symmetry breaking may only occur on one side, a map that appears in one
may not take this symmetry breaking into account.
:type r_matcher: GraphMatcher instance
:param r_matcher: iterator over isomorphic mappings of the input reactant to
template reactant
:type p_matcher: GraphMatcher instance
:param p_matcher: iterator over isomorphic mappings of the input product to
template product
:type input_rc: list of ints
:param input_rc: atom indices of input reaction center
:type template_rc: list of ints
:param template_rc: atom indices of template reaction center
"""
def image(mapping, domain):
"""
return image of a map given domain
"""
return set(v for k, v in mapping.items() if k in domain)
p_maps = set(
frozenset(isomorph.items())
for isomorph in p_matcher.isomorphisms_iter()
if image(isomorph, input_rc) == set(template_rc))
for isomorph in r_matcher.isomorphisms_iter():
if frozenset(isomorph.items()) in p_maps:
yield isomorph