com.actelion.research.chem.reaction.mapping.RootAtomPairSource Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
/*
* Copyright (c) 1997 - 2016
* Actelion Pharmaceuticals Ltd.
* Gewerbestrasse 16
* CH-4123 Allschwil, Switzerland
*
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of the the copyright holder nor the
* names of its contributors may be used to endorse or promote products
* derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* @author Thomas Sander
*/
package com.actelion.research.chem.reaction.mapping;
import com.actelion.research.chem.Canonizer;
import com.actelion.research.chem.CanonizerBaseValue;
import com.actelion.research.chem.Molecule;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.util.ByteArrayComparator;
import com.actelion.research.util.IntArrayComparator;
import java.nio.charset.StandardCharsets;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.TreeMap;
import java.util.TreeSet;
public class RootAtomPairSource {
private static final int MAX_PAIR_SEQUENCES = 64;
private static final int MAX_ENVIRONMENT_RADIUS = 7;
private static final int MIN_ENVIRONMENT_RADIUS = 2;
private static final int PSEUDO_MAP_NO_SKIPPED_PAIR = -99999;
private ArrayList mPairBuffer;
private StereoMolecule mReactant,mProduct;
private Canonizer mReactantCanonizer,mProductCanonizer;
private CanonizerBaseValue[] mCanBase;
private int mSequenceCount,mCurrentRadius,mManualMapCount,mMappableAtomCount,mCurrentEnvIndex0,mCurrentEnvIndex1,mCurrentEnvIndex2,mCurrentEnvIndex3;
private RootAtomPairDecisionHelper mDecisionHelper;
private boolean mIsStoichiometric;
private int[] mReactantRank,mProductRank,mReactantFragmentNo,mProductFragmentNo,mReactantMapNo,mProductMapNo;
private int mAtomBits,mMaxConnAtoms,mHighestReactionRank,mHighestProductRank;
private boolean[] mReactantFragmentUsed,mProductFragmentUsed;
private TreeMap[] mEnvToAtomsMap; // index: radius
private byte[][][] mEnvKey;
public RootAtomPairSource(StereoMolecule reactant, StereoMolecule product, int[] reactantMapNo, int[] productMapNo) {
mReactant = reactant;
mProduct = product;
mReactantMapNo = reactantMapNo;
mProductMapNo = productMapNo;
// Reactant and product atom ranks are updated with every assigned pair in order to reflect
// decreasing symmetry due to unsymmetrically assigned atom mapping.
initializeAtomRanking();
initializeRanks();
mEnvToAtomsMap = buildEnvToAtomMaps();
mEnvKey = getEnvKeys(mEnvToAtomsMap);
mMappableAtomCount = 0;
for (byte[] envKey:mEnvToAtomsMap[0].keySet()) {
int[][] atoms = mEnvToAtomsMap[0].get(envKey);
mMappableAtomCount += Math.min(atoms[0].length, atoms[1].length);
}
mIsStoichiometric = mMappableAtomCount == mReactant.getAtoms()
&& mReactant.getAtoms() == mProduct.getAtoms();
mReactantFragmentNo = new int[mReactant.getAtoms()];
mProductFragmentNo = new int[mProduct.getAtoms()];
mReactantFragmentUsed = new boolean[mReactant.getFragmentNumbers(mReactantFragmentNo, false, false)];
mProductFragmentUsed = new boolean[mProduct.getFragmentNumbers(mProductFragmentNo, false, false)];
mManualMapCount = assignManuallyMappedPairs();
mCurrentRadius = MAX_ENVIRONMENT_RADIUS;
mSequenceCount = 0;
}
private void initializeAtomRanking() {
int maxAtomCount = Math.max(mReactant.getAtoms(), mProduct.getAtoms());
mAtomBits = Canonizer.getNeededBits(maxAtomCount);
mMaxConnAtoms = 2;
for (int atom=0; atom[] buildEnvToAtomMaps() {
byte[][][] reactantAtomEnv = classifyAtomEnvironment(mReactant);
byte[][][] productAtomEnv = classifyAtomEnvironment(mProduct);
TreeMap[] envMap = new TreeMap[MAX_ENVIRONMENT_RADIUS+1];
TreeMap[] reactantMap = buildEnvToAtomMaps(mReactant, reactantAtomEnv);
TreeMap[] productMap = buildEnvToAtomMaps(mProduct, productAtomEnv);
for (int radius=0; radius<=MAX_ENVIRONMENT_RADIUS; radius++) {
envMap[radius] = new TreeMap<>(new ByteArrayComparator());
for (byte[] envKey:reactantMap[radius].keySet()) {
int[] reactantAtoms = reactantMap[radius].get(envKey);
int[] productAtoms = productMap[radius].get(envKey);
if (productAtoms != null) {
int[][] atoms = new int[2][];
atoms[0] = reactantAtoms;
atoms[1] = productAtoms;
envMap[radius].put(envKey, atoms);
}
}
}
return envMap;
}
private TreeMap[] buildEnvToAtomMaps(StereoMolecule mol, byte[][][] atomEnv) {
TreeMap[] map = new TreeMap[MAX_ENVIRONMENT_RADIUS+1];
for (int radius=0; radius<=MAX_ENVIRONMENT_RADIUS; radius++) {
map[radius] = new TreeMap<>(new ByteArrayComparator());
for (int atom=0; atom[] envMap) {
byte[][][] keys = new byte[MAX_ENVIRONMENT_RADIUS+1][][];
for (int radius=0; radius<=MAX_ENVIRONMENT_RADIUS; radius++) {
keys[radius] = new byte[envMap[radius].size()][];
int index = 0;
for (byte[] key:envMap[radius].keySet())
keys[radius][index++] = key;
}
return keys;
}
/**
* Assigns first mapping numbers to manually mapped atoms and add the pair list to the buffer.
* @return list of manually mapped atom pairs
*/
private int assignManuallyMappedPairs() {
mPairBuffer = new ArrayList<>();
int mapNo = 1;
int maxMapNo = 0;
for (int atom=0; atom
* - they match in terms of circular fragment on reactant and products side
* - if multiple symmetrically equivalent pairs exist, exactly one of them is marked as allowed root pair
* - each of the reactant and product atoms have at least one unmapped neighbour
* @return pair of currently obvious root atom pairs
*/
public RootAtomPair nextPair() {
RootAtomPair pair = nextRawPair();
while (pair != null) {
boolean reactantAtomHasUnmappedNeighbours = false;
for (int i=0; i= 0) {
// We create starting pairs of reasonably similar atoms with similarity derived priority
while (mCurrentRadius >= MIN_ENVIRONMENT_RADIUS
&& mCurrentEnvIndex0 < mEnvKey[mCurrentRadius].length) {
byte[] envKey = mEnvKey[mCurrentRadius][mCurrentEnvIndex0];
int[][] atoms = mEnvToAtomsMap[mCurrentRadius].get(envKey);
if (mReactant.getAtomicNo(atoms[0][0]) == 6) {
RootAtomPair pair = makePairsFromSimilarAtoms(atoms[0], atoms[1]);
if (pair != null)
return pair;
}
mCurrentEnvIndex0++; // go to next environment key once all potential pairs with current environment are depleted
}
while (mCurrentRadius >= MIN_ENVIRONMENT_RADIUS
&& mCurrentEnvIndex1 < mEnvKey[mCurrentRadius].length) {
byte[] envKey = mEnvKey[mCurrentRadius][mCurrentEnvIndex1];
int[][] atoms = mEnvToAtomsMap[mCurrentRadius].get(envKey);
if (mReactant.getAtomicNo(atoms[0][0]) != 6) {
// with equal environment size, we prefer carbon root atoms TODO
RootAtomPair pair = makePairsFromSimilarAtoms(atoms[0], atoms[1]);
if (pair != null)
return pair;
}
mCurrentEnvIndex1++; // go to next environment key once all potential pairs with current environment are depleted
}
// We create a low priority starting pair if we just have one atom of a kind
while (mIsStoichiometric && mCurrentRadius == 0
&& mCurrentEnvIndex2 < mEnvKey[0].length) {
byte[] envKey = mEnvKey[0][mCurrentEnvIndex2++];
int[][] atoms = mEnvToAtomsMap[mCurrentRadius].get(envKey);
if (atoms[0].length == 1 && atoms[1].length == 1) {
RootAtomPair pair = tryCreatePair(atoms[0][0], atoms[1][0]);
if (pair != null)
return pair;
}
}
// Check with decreasing sphere size down to atomicNo level, whether all reactant or product atoms
// are in separate fragments and are equivalent. If this is the case, and if we have at least
// as many equivalent atoms as atoms on the other side, then we just match reactant atoms
// to product atoms in order of appearance.
while (mCurrentEnvIndex3 < mEnvKey[mCurrentRadius].length) {
byte[] envKey = mEnvKey[mCurrentRadius][mCurrentEnvIndex3++];
int[][] atoms = mEnvToAtomsMap[mCurrentRadius].get(envKey);
if (atoms[0].length == 1 && atoms[1].length == 1) {
RootAtomPair pair = tryCreatePair(atoms[0][0], atoms[1][0]);
if (pair != null)
return pair;
}
else if ((atoms[0].length >= atoms[1].length
&& areInDistinctEquivalentFragments(atoms[0], mReactantFragmentNo, mReactantFragmentUsed, mReactantRank))
|| (atoms[1].length >= atoms[0].length
&& areInDistinctEquivalentFragments(atoms[1], mProductFragmentNo, mProductFragmentUsed, mProductRank))) {
for (int i=0; i rankPairSet = new TreeSet<>(new IntArrayComparator());
for (int ra:reactantAtom) {
for (int pa:productAtom) {
int[] pair = new int[2];
pair[0] = mReactantRank[ra];
pair[1] = mProductRank[pa];
rankPairSet.add(pair);
}
}
return rankPairSet.toArray(new int[0][]);
}
/**
* Increases the symmetry rank of the specified atom in regard to all other atoms with
* the same rank. This split of formerly equal ranks accounts for the asymmetry introduced
* by matching this atom to another atom on the other side of the reaction.
* The updated rank is then propagated recursively through all neighbours.
* This, of course, happens only, if the atom's initial rank is shared by at least another atom.
* @param mol
* @param atomRank reactant's or product's atom ranks to be modified in place
* @param atom
* @return
*/
private int elevateAtomRank(StereoMolecule mol, int[] atomRank, int atom, int maxRank) {
int currentAtomRank = atomRank[atom];
boolean sharedRankFound = false;
for (int i=0; i currentAtomRank)
atomRank[i]++;
int oldMaxRank;
maxRank++;
do {
oldMaxRank = maxRank;
canCalcNextBaseValues(mol, atomRank);
maxRank = consolidateRanking(atomRank);
} while (oldMaxRank != maxRank);
return maxRank;
}
private void canCalcNextBaseValues(StereoMolecule mol, int[] atomRank) {
int[] connRank = new int[mMaxConnAtoms];
for (int atom=0; atom=mol.getAllConnAtoms(atom)) {
int rank = 2 * atomRank[mol.getConnAtom(atom, i)];
int connBond = mol.getConnBond(atom, i);
if (mol.getBondOrder(connBond) == 2)
if (!mol.isAromaticBond(connBond))
rank++; // set a flag for non-aromatic double bond
int j;
for (j = 0; j < neighbour; j++)
if (rank < connRank[j])
break;
for (int k = neighbour; k > j; k--)
connRank[k] = connRank[k - 1];
connRank[j] = rank;
neighbour++;
}
}
mCanBase[atom].init(atom);
mCanBase[atom].add(mAtomBits, atomRank[atom]);
for (int i=neighbours; i {
public int reactantAtom,productAtom;
public RootAtomPair(int reactantAtom, int productAtom) {
this.reactantAtom = reactantAtom;
this.productAtom = productAtom;
}
@Override
public int compareTo(RootAtomPair pair) {
return this.reactantAtom < pair.reactantAtom ? -1
: this.reactantAtom > pair.reactantAtom ? 1
: this.productAtom < pair.productAtom ? -1
: this.productAtom > pair.productAtom ? 1 : 0;
}
}