feat: add getHoseFragment on TopicMolecule

package.json

@@ -40,23 +40,23 @@
   },
   "homepage": "https://github.com/cheminfo/openchemlib-utils#readme",
   "devDependencies": {
-    "@types/node": "^20.12.6",
-    "@vitest/coverage-v8": "^1.4.0",
+    "@types/node": "^20.12.11",
+    "@vitest/coverage-v8": "^1.6.0",
     "cheminfo-build": "^1.2.0",
     "cheminfo-types": "^1.7.3",
     "eslint": "^8.57.0",
-    "eslint-config-cheminfo-typescript": "^12.2.0",
+    "eslint-config-cheminfo-typescript": "^12.4.0",
     "fifo-logger": "^1.0.0",
-    "mf-parser": "^3.1.0",
+    "mf-parser": "^3.1.1",
     "openchemlib": "^8.9.0",
     "prettier": "^3.2.5",
-    "rimraf": "^5.0.5",
+    "rimraf": "^5.0.7",
     "ts-node": "^10.9.2",
-    "typescript": "^5.4.4",
-    "vitest": "^1.4.0"
+    "typescript": "^5.4.5",
+    "vitest": "^1.6.0"
   },
   "dependencies": {
-    "atom-sorter": "^2.0.0",
+    "atom-sorter": "^2.0.1",
     "ensure-string": "^1.2.0",
     "get-value": "^3.0.1",
     "ml-floyd-warshall": "^3.0.1",

src/hose/getHoseCodesForAtomsInternal.js

@@ -39,8 +39,8 @@ export function getHoseCodesForAtomsInternal(molecule, options = {}) {
   const results = [];
   let min = 0;
   let max = 0;
-  const atomMask = new Array(molecule.getAllAtoms());
-  const atomList = new Array(molecule.getAllAtoms());
+  const atomMask = new Uint8Array(molecule.getAllAtoms());
+  const atomList = new Uint8Array(molecule.getAllAtoms());
 
   for (let sphere = 0; sphere <= maxSphereSize; sphere++) {
     if (max === 0) {

src/path/getAllAtomsPaths.ts

@@ -8,7 +8,7 @@ export interface GetAllAtomsPathsOptions {
   maxPathLength?: number;
 }
 
-interface AtomPath {
+export interface AtomPath {
   path: number[];
   distance: number;
 }

src/topic/TopicMolecule.ts

@@ -1,6 +1,8 @@
 import type { Molecule } from 'openchemlib';
 
+import { AtomPath, getAllAtomsPaths } from '../path/getAllAtomsPaths.js';
 import { getConnectivityMatrix } from '../util/getConnectivityMatrix.js';
+import { tagAtom } from '../util/tagAtom.js';
 
 import { HoseCodesOptions } from './HoseCodesOptions.js';
 import { getCanonizedDiaIDs } from './getCanonizedDiaIDs';
@@ -17,6 +19,34 @@ interface ToMolfileOptions {
   version?: 2 | 3;
 }
 
+interface TopicMoleculeOptions extends HoseCodesOptions {
+  maxPathLength?: 5;
+}
+
+interface GetAtomPathOptions {
+  /*
+   * The distance between the two atoms. If not specified, all the distances will be considered
+   */
+  distance?: number;
+}
+
+interface GetHoseFragmentOptions {
+  /*
+   * The sphere size around any selected atoms to consider. Default is 2
+   */
+  sphereSize?: number;
+  /**
+   * The atoms to tag in the fragment
+   * @default rootAtoms
+   */
+  tagAtoms?: number[];
+  /**
+   * The function to tag the atoms in place !
+   * @default tagAtom
+   */
+  tagAtomFct?: (molecule: Molecule, iAtom: number) => undefined;
+}
+
 /**
  * This class deals with topicity information and hose codes
  * It is optimized to avoid recalculation of the same information
@@ -25,13 +55,13 @@ export class TopicMolecule {
   private readonly originalMolecule: Molecule;
   molecule: Molecule;
   idCode: string;
-  options: HoseCodesOptions;
+  options: TopicMoleculeOptions;
 
   private cache: any;
 
-  constructor(molecule: Molecule, options: HoseCodesOptions = {}) {
+  constructor(molecule: Molecule, options: TopicMoleculeOptions = {}) {
     this.originalMolecule = molecule;
-    this.options = options;
+    this.options = { maxPathLength: 5, ...options };
     this.idCode = molecule.getIDCode();
     this.molecule = this.originalMolecule.getCompactCopy();
     this.molecule.ensureHelperArrays(
@@ -69,7 +99,99 @@ export class TopicMolecule {
     }
   }
 
-  toMolfile(options: ToMolfileOptions = {}) {
+  getHoseFragment(
+    rootAtoms: number[],
+    options: GetHoseFragmentOptions = {},
+  ): Molecule {
+    const {
+      sphereSize = 2,
+      tagAtoms = rootAtoms,
+      tagAtomFct = tagAtom,
+    } = options;
+    this.moleculeWithH.ensureHelperArrays(
+      this.moleculeWithH.getOCL().Molecule.cHelperNeighbours,
+    );
+
+    const copy = this.moleculeWithH.getCompactCopy();
+    copy.ensureHelperArrays(copy.getOCL().Molecule.cHelperNeighbours);
+
+    for (let i = 0; i < copy.getAllAtoms(); i++) {
+      copy.setAtomMass(i, copy.getAtomMass(i) + 2);
+    }
+
+    const atomMask = new Array(copy.getAllAtoms()).fill(false);
+    const atomList = new Uint8Array(copy.getAllAtoms());
+    const atomMapping = new Array(copy.getAllAtoms()).fill(-1);
+    const Molecule = copy.getOCL().Molecule;
+    const fragment = new Molecule(0, 0);
+    fragment.setFragment(true);
+    let min = 0;
+    let max = 0;
+    for (let sphere = 0; sphere <= sphereSize; sphere++) {
+      if (max === 0) {
+        for (const rootAtom of rootAtoms) {
+          atomList[max] = rootAtom;
+          atomMask[rootAtom] = 1;
+          max++;
+        }
+      } else {
+        let newMax = max;
+        for (let i = min; i < max; i++) {
+          const atom = atomList[i];
+          for (let j = 0; j < this.moleculeWithH.getAllConnAtoms(atom); j++) {
+            const connAtom = this.moleculeWithH.getConnAtom(atom, j);
+            if (!atomMask[connAtom]) {
+              atomMask[connAtom] = 1;
+              atomList[newMax++] = connAtom;
+            }
+          }
+        }
+        min = max;
+        max = newMax;
+      }
+    }
+    copy.copyMoleculeByAtoms(fragment, atomMask, true, atomMapping);
+    for (let i = 0; i < fragment.getAllAtoms(); i++) {
+      fragment.setAtomMass(i, fragment.getAtomMass(i) - 2);
+    }
+
+    for (const atom of tagAtoms) {
+      tagAtomFct(fragment, atomMapping[atom]);
+    }
+
+    return fragment;
+  }
+
+  getAtomPaths(atom1: number, atom2: number, options: GetAtomPathOptions = {}) {
+    const { distance } = options;
+    if (distance !== undefined && distance > this.options.maxPathLength) {
+      throw new Error(
+        'The distance is too long, you should increase the maxPathLength when instanciating the TopicMolecule',
+      );
+    }
+    const atomPaths = this.atomsPaths[atom1];
+    const minDistance = distance || 0;
+    const maxDistance = distance || this.options.maxPathLength;
+    const paths = [];
+    for (let i = minDistance; i <= maxDistance; i++) {
+      for (const atomPath of atomPaths[i]) {
+        if (atomPath.path.at(-1) === atom2) {
+          paths.push(atomPath.path);
+        }
+      }
+    }
+    return paths;
+  }
+
+  get atomsPaths(): AtomPath[][][] {
+    if (this.cache.atomsPaths) return this.cache.atomsPaths;
+    this.cache.atomsPaths = getAllAtomsPaths(this.moleculeWithH, {
+      maxPathLength: this.options.maxPathLength,
+    });
+    return this.cache.atomsPaths;
+  }
+
+  toMolfile(options: ToMolfileOptions = {}): string {
     const { version = 2 } = options;
     if (version === 2) {
       return this.molecule.toMolfile();
@@ -105,13 +227,13 @@ export class TopicMolecule {
   /**
    * Returns a molecule with all the hydrogens added. The order is NOT canonized
    */
-  get moleculeWithH() {
+  get moleculeWithH(): Molecule {
     if (this.cache.moleculeWithH) return this.cache.moleculeWithH;
     this.cache.moleculeWithH = getMoleculeWithH(this.molecule);
     return this.cache.moleculeWithH;
   }
 
-  private get xMolecule() {
+  private get xMolecule(): Molecule {
     if (this.cache.xMolecule) return this.cache.xMolecule;
     this.cache.xMolecule = getXMolecule(this.moleculeWithH);
     return this.cache.xMolecule;

src/topic/__tests__/TopicMolecule.pathHose.test.js

@@ -0,0 +1,40 @@
+import { Molecule } from "openchemlib";
+import { test, expect } from "vitest";
+
+import { makeRacemic } from "../../util/makeRacemic.js";
+import { TopicMolecule } from "../TopicMolecule.js";
+
+
+test("TopicMolecule.path", () => {
+
+  const molecule = Molecule.fromSmiles('CCCCCO');
+  const topicMolecule = new TopicMolecule(molecule);
+
+  expect(topicMolecule.getAtomPaths(0, 2, { distance: 2 })).toStrictEqual([[0, 1, 2]])
+  expect(topicMolecule.getAtomPaths(0, 2, { distance: 3 })).toStrictEqual([])
+
+})
+
+test("TopicMolecule.getHoseFragment", async () => {
+  const molecule = Molecule.fromSmiles('Cl/C=C/CCCCl');
+  //const molecule = Molecule.fromSmiles('C1=CC=CC=C1');
+  molecule.addImplicitHydrogens();
+
+  const topicMolecule = new TopicMolecule(molecule);
+
+  const fragment = topicMolecule.getHoseFragment([0, 1, 2, 9], {
+    sphereSize: 1, tagAtoms: [2, 9],
+    tagAtomFct: (fragment, atom) => {
+      fragment.setAtomMass(atom, 255)
+    }
+  })
+
+  makeRacemic(fragment);
+
+
+
+  expect(fragment.toIsomericSmiles().replaceAll('255', '*')).toBe("[*H]C/[*C]=C/Cl")
+  expect(fragment.toIsomericSmiles({ createSmarts: true }).replaceAll('255', '*')).toBe('[*H]C/[*C;!H0]=[C;!H0]/Cl')
+  expect(fragment.toIsomericSmiles({ kekulizedOutput: true }).replaceAll('255', '*')).toBe('[*H]C/[*C]=C/Cl')
+
+})