From c213a7d3c3b3b7c0c0ee2724d7d9ea42326fbf6d Mon Sep 17 00:00:00 2001 From: Jose Duarte Date: Tue, 1 Jan 2019 16:39:45 -0800 Subject: [PATCH 01/11] Neighbor indices now uses spatial hashing --- .../nbio/structure/asa/AsaCalculator.java | 71 +++++++++++++++++-- .../nbio/structure/asa/TestAsaCalc.java | 65 +++++++++++++++-- 2 files changed, 127 insertions(+), 9 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index dc34ff8782..224221c4f8 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -21,12 +21,13 @@ package org.biojava.nbio.structure.asa; import org.biojava.nbio.structure.*; +import org.biojava.nbio.structure.contact.Contact; +import org.biojava.nbio.structure.contact.Grid; import org.slf4j.Logger; import org.slf4j.LoggerFactory; import javax.vecmath.Point3d; -import java.util.ArrayList; -import java.util.TreeMap; +import java.util.*; import java.util.concurrent.ExecutorService; import java.util.concurrent.Executors; @@ -101,6 +102,7 @@ public void run() { private int nThreads; private Point3d[] spherePoints; private double cons; + private List contacts; /** * Constructs a new AsaCalculator. Subsequently call {@link #calculateAsas()} @@ -319,8 +321,9 @@ private Point3d[] generateSpherePoints(int nSpherePoints) { /** * Returns list of indices of atoms within probe distance to atom k. * @param k index of atom for which we want neighbor indices + * @return the indices of neighboring atoms */ - private Integer[] findNeighborIndices(int k) { + Integer[] findNeighborIndices(int k) { // looking at a typical protein case, number of neighbours are from ~10 to ~50, with an average of ~30 // Thus 40 seems to be a good compromise for the starting capacity ArrayList neighbor_indices = new ArrayList<>(40); @@ -343,9 +346,69 @@ private Integer[] findNeighborIndices(int k) { return indicesArray; } + /** + * Returns list of indices of atoms within probe distance to atom k, + * using spatial hashing to avoid all to all distance calculation. + * @param k index of atom for which we want neighbor indices + * @return the indices of neighboring atoms + */ + Integer[] findNeighborIndicesSpatialHashing(int k) { + + if (contacts == null) { + contacts = calcContacts(); + } + + // looking at a typical protein case, number of neighbours are from ~10 to ~50, with an average of ~30 + // Thus 40 seems to be a good compromise for the starting capacity + ArrayList neighbor_indices = new ArrayList<>(40); + + double radius = radii[k] + probe + probe; + + for (Contact contact : contacts) { + double dist = contact.getDistance(); + int i; + if (contact.getJ() == k) { + i = contact.getI(); + } else if (contact.getI() == k) { + i = contact.getJ(); + } else { + continue; + } + if (dist < radius + radii[i]) { + neighbor_indices.add(i); + } + } + + Integer[] indicesArray = new Integer[neighbor_indices.size()]; + indicesArray = neighbor_indices.toArray(indicesArray); + return indicesArray; + } + + Point3d[] getAtomCoords() { + return atomCoords; + } + + private List calcContacts() { + double maxRadius = maxValue(radii); + double cutoff = maxRadius + maxRadius + probe + probe; + Grid grid = new Grid(cutoff); + grid.addCoords(atomCoords); + return grid.getIndicesContacts(); + } + + private static double maxValue(double[] array) { + double max = array[0]; + for (int i = 0; i < array.length; i++) { + if (array[i] > max) { + max = array[i]; + } + } + return max; + } + private double calcSingleAsa(int i) { Point3d atom_i = atomCoords[i]; - Integer[] neighbor_indices = findNeighborIndices(i); + Integer[] neighbor_indices = findNeighborIndicesSpatialHashing(i); int n_neighbor = neighbor_indices.length; int j_closest_neighbor = 0; double radius = probe + radii[i]; diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java index e7ea1d2481..ea163fa212 100644 --- a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java +++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java @@ -20,16 +20,17 @@ */ package org.biojava.nbio.structure.asa; -import junit.framework.TestCase; import org.biojava.nbio.structure.Structure; import org.biojava.nbio.structure.StructureException; import org.biojava.nbio.structure.StructureIO; import org.biojava.nbio.structure.io.mmcif.ChemCompGroupFactory; import org.biojava.nbio.structure.io.mmcif.DownloadChemCompProvider; -import org.junit.Assert; +import static org.junit.Assert.*; import org.junit.Test; import java.io.IOException; +import java.util.ArrayList; +import java.util.List; /** * Testing of Accessible Surface Area calculations @@ -70,12 +71,66 @@ public void testAsa3PIU() throws StructureException, IOException { //System.out.println(groupAsa.getGroup().getPDBName() + " " + groupAsa.getGroup().getResidueNumber() + " " + groupAsa.getAsaU()); totResidues+=groupAsa.getAsaU(); - Assert.assertTrue(groupAsa.getRelativeAsaU() <= 1.0); + assertTrue(groupAsa.getRelativeAsaU() <= 1.0); } - Assert.assertEquals(totAtoms, totResidues, 0.000001); + assertEquals(totAtoms, totResidues, 0.000001); - Assert.assertEquals(17462.0, totAtoms, 1.0); + assertEquals(17462.0, totAtoms, 1.0); + + } + + @Test + public void testNeighborIndicesFinding() throws StructureException, IOException { + // important: without this the tests can fail when running in maven (but not in IDE) + // that's because it depends on the order on how tests were run - JD 2018-03-10 + ChemCompGroupFactory.setChemCompProvider(new DownloadChemCompProvider()); + + Structure structure = StructureIO.getStructure("3PIU"); + + AsaCalculator asaCalc = new AsaCalculator(structure, + AsaCalculator.DEFAULT_PROBE_SIZE, + 1000, 1, false); + + for (int indexToTest =0; indexToTest < asaCalc.getAtomCoords().length; indexToTest++) { + //int indexToTest = 198; + + Integer[] nbsSh = asaCalc.findNeighborIndicesSpatialHashing(indexToTest); + + Integer[] nbs = asaCalc.findNeighborIndices(indexToTest); + + int countNotInNbs = 0; + List listOfMatchingIndices = new ArrayList<>(); + for (int i = 0; i < nbsSh.length; i++) { + boolean contained = false; + for (int j = 0; j < nbs.length; j++) { + if (nbs[j].equals(nbsSh[i])) { + listOfMatchingIndices.add(j); + contained = true; + break; + } + } + if (!contained) { + countNotInNbs++; + } + } + + //System.out.println("In nbsSh but not in nbs: " + countNotInNbs); + //System.out.println("Number of matching indices: " + listOfMatchingIndices.size()); + +// for (int i = 0; i Date: Tue, 1 Jan 2019 16:56:45 -0800 Subject: [PATCH 02/11] Supporting neighbors with and without SH --- .../nbio/structure/asa/AsaCalculator.java | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 224221c4f8..0d75c2812f 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -104,6 +104,8 @@ public void run() { private double cons; private List contacts; + private boolean useSpatialHashingForNeighbors; + /** * Constructs a new AsaCalculator. Subsequently call {@link #calculateAsas()} * or {@link #getGroupAsas()} to calculate the ASAs @@ -122,6 +124,8 @@ public AsaCalculator(Structure structure, double probe, int nSpherePoints, int n this.probe = probe; this.nThreads = nThreads; + this.useSpatialHashingForNeighbors = true; + // initialising the radii by looking them up through AtomRadii radii = new double[atomCoords.length]; for (int i=0;i Date: Tue, 1 Jan 2019 17:47:46 -0800 Subject: [PATCH 03/11] A test, demoing that SH actually performs much worse --- .../nbio/structure/asa/AsaCalculator.java | 8 ++- .../nbio/structure/asa/TestAsaCalc.java | 61 +++++++++++++++++-- 2 files changed, 62 insertions(+), 7 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 0d75c2812f..215c91a3bd 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -59,6 +59,8 @@ public class AsaCalculator { public static final double DEFAULT_PROBE_SIZE = 1.4; public static final int DEFAULT_NTHREADS = 1; + public static final boolean DEFAULT_USE_SPATIAL_HASHING = true; + // Chothia's amino acid atoms vdw radii @@ -124,7 +126,7 @@ public AsaCalculator(Structure structure, double probe, int nSpherePoints, int n this.probe = probe; this.nThreads = nThreads; - this.useSpatialHashingForNeighbors = true; + this.useSpatialHashingForNeighbors = DEFAULT_USE_SPATIAL_HASHING; // initialising the radii by looking them up through AtomRadii radii = new double[atomCoords.length]; @@ -154,7 +156,7 @@ public AsaCalculator(Atom[] atoms, double probe, int nSpherePoints, int nThreads this.probe = probe; this.nThreads = nThreads; - this.useSpatialHashingForNeighbors = true; + this.useSpatialHashingForNeighbors = DEFAULT_USE_SPATIAL_HASHING; for (Atom atom:atoms) { if (atom.getElement()==Element.H) @@ -196,7 +198,7 @@ public AsaCalculator(Point3d[] atomCoords, double probe, int nSpherePoints, int this.probe = probe; this.nThreads = nThreads; - this.useSpatialHashingForNeighbors = true; + this.useSpatialHashingForNeighbors = DEFAULT_USE_SPATIAL_HASHING; // initialising the radii to the given radius for all atoms radii = new double[atomCoords.length]; diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java index ea163fa212..e4cf36c7e2 100644 --- a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java +++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java @@ -20,12 +20,11 @@ */ package org.biojava.nbio.structure.asa; -import org.biojava.nbio.structure.Structure; -import org.biojava.nbio.structure.StructureException; -import org.biojava.nbio.structure.StructureIO; +import org.biojava.nbio.structure.*; import org.biojava.nbio.structure.io.mmcif.ChemCompGroupFactory; import org.biojava.nbio.structure.io.mmcif.DownloadChemCompProvider; import static org.junit.Assert.*; + import org.junit.Test; import java.io.IOException; @@ -36,7 +35,7 @@ * Testing of Accessible Surface Area calculations * * - * @author duarte_j + * @author Jose Duarte * */ public class TestAsaCalc { @@ -133,4 +132,58 @@ public void testNeighborIndicesFinding() throws StructureException, IOException } } + + @Test + public void testPerformance() throws StructureException, IOException { + // important: without this the tests can fail when running in maven (but not in IDE) + // that's because it depends on the order on how tests were run - JD 2018-03-10 + ChemCompGroupFactory.setChemCompProvider(new DownloadChemCompProvider()); + + Structure structure = StructureIO.getStructure("4F5X"); + Chain c = structure.getPolyChainByPDB("W"); + Atom[] atoms = StructureTools.getAllAtomArray(c); + System.out.printf("Total of %d atoms\n", atoms.length); + + int nThreads = 1; + // 1. WITH SPATIAL HASHING + + long start = System.currentTimeMillis(); + AsaCalculator asaCalc = new AsaCalculator(atoms, + AsaCalculator.DEFAULT_PROBE_SIZE, + 100, nThreads); + asaCalc.setUseSpatialHashingForNeighbors(true); + + double[] asas = asaCalc.calculateAsas(); + long end = System.currentTimeMillis(); + System.out.printf("ASA calculation took %6.2f s with spatial hashing\n", (end-start)/1000.0); + + double totAtoms = 0; + for (double asa:asas) { + totAtoms += asa; + } + double withSH = totAtoms; + System.out.printf("Total ASA is %6.2f \n", totAtoms); + + + // 2. WITHOUT SPATIAL HASHING + start = System.currentTimeMillis(); + asaCalc = new AsaCalculator(atoms, + AsaCalculator.DEFAULT_PROBE_SIZE, + 100, nThreads); + asaCalc.setUseSpatialHashingForNeighbors(false); + + asas = asaCalc.calculateAsas(); + end = System.currentTimeMillis(); + System.out.printf("ASA calculation took %6.2f s without spatial hashing\n", (end-start)/1000.0); + + totAtoms = 0; + for (double asa:asas) { + totAtoms += asa; + } + double withoutSH = totAtoms; + System.out.printf("Total ASA is %6.2f \n", totAtoms); + + assertEquals(withoutSH, withSH, 0.000001); + + } } From 11002993661bdab0410dbc7a6d651e0ae7f4aea8 Mon Sep 17 00:00:00 2001 From: Jose Duarte Date: Tue, 1 Jan 2019 18:56:02 -0800 Subject: [PATCH 04/11] Doing all neighbor indices upfront --- .../nbio/structure/asa/AsaCalculator.java | 111 ++++++++++-------- .../nbio/structure/asa/TestAsaCalc.java | 18 ++- 2 files changed, 71 insertions(+), 58 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 215c91a3bd..22b28a47e6 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -104,7 +104,7 @@ public void run() { private int nThreads; private Point3d[] spherePoints; private double cons; - private List contacts; + private int[][] neighborIndices; private boolean useSpatialHashingForNeighbors; @@ -250,6 +250,12 @@ public double[] calculateAsas() { double[] asas = new double[atomCoords.length]; + if (useSpatialHashingForNeighbors) { + neighborIndices = findNeighborIndicesSpatialHashing(); + } else { + neighborIndices = findNeighborIndices(); + } + if (nThreads<=1) { // (i.e. it will also be 1 thread if 0 or negative number specified) for (int i=0;i neighbor_indices = new ArrayList<>(40); + int[][] findNeighborIndices() { - double radius = radii[k] + probe + probe; + int[][] nbsIndices = new int[atomCoords.length][]; - for (int i=0;i thisNbIndices = new ArrayList<>(); - if (dist < radius + radii[i]) { - neighbor_indices.add(i); + for (int i = 0; i < atomCoords.length; i++) { + if (i == k) continue; + + double dist = atomCoords[i].distance(atomCoords[k]); + + if (dist < radius + radii[i]) { + thisNbIndices.add(i); + } } + int[] indicesArray = new int[thisNbIndices.size()]; + for (int i=0;i neighbor_indices = new ArrayList<>(40); + List contactList = calcContacts(); - double radius = radii[k] + probe + probe; + for (int k=0; k thisNbIndices = new ArrayList<>(); + + // TODO make this the outer loop + for (Contact contact : contactList) { + double dist = contact.getDistance(); + int i; + if (contact.getJ() == k) { + i = contact.getI(); + } else if (contact.getI() == k) { + i = contact.getJ(); + } else { + continue; + } + if (dist < radius + radii[i]) { + thisNbIndices.add(i); + } } + + int[] indicesArray = new int[thisNbIndices.size()]; + for (int i=0;i calcContacts() { double maxRadius = maxValue(radii); double cutoff = maxRadius + maxRadius + probe + probe; + logger.debug("Max radius is {}, cutoff is {}", maxRadius, cutoff); Grid grid = new Grid(cutoff); grid.addCoords(atomCoords); return grid.getIndicesContacts(); @@ -422,13 +433,9 @@ private static double maxValue(double[] array) { private double calcSingleAsa(int i) { Point3d atom_i = atomCoords[i]; - Integer[] neighbor_indices; - if (useSpatialHashingForNeighbors) { - neighbor_indices = findNeighborIndicesSpatialHashing(i); - } else { - neighbor_indices = findNeighborIndices(i); - } - int n_neighbor = neighbor_indices.length; + + int n_neighbor = neighborIndices[i].length; + int[] neighbor_indices = neighborIndices[i]; int j_closest_neighbor = 0; double radius = probe + radii[i]; @@ -579,7 +586,7 @@ private static double getRadiusForNucl(NucleotideImpl nuc, Atom atom) { * * If atom is neither part of a nucleotide nor of a standard aminoacid, * the default vdw radius for the element is returned. If atom is of - * unknown type (element) the vdw radius of {@link #Element().N} is returned + * unknown type (element) the vdw radius of {@link Element().N} is returned * * @param atom * @return diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java index e4cf36c7e2..1f6e4696b4 100644 --- a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java +++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java @@ -91,19 +91,21 @@ public void testNeighborIndicesFinding() throws StructureException, IOException AsaCalculator.DEFAULT_PROBE_SIZE, 1000, 1, false); - for (int indexToTest =0; indexToTest < asaCalc.getAtomCoords().length; indexToTest++) { - //int indexToTest = 198; + int[][] allNbsSh = asaCalc.findNeighborIndicesSpatialHashing(); - Integer[] nbsSh = asaCalc.findNeighborIndicesSpatialHashing(indexToTest); + int[][] allNbs = asaCalc.findNeighborIndices(); - Integer[] nbs = asaCalc.findNeighborIndices(indexToTest); + for (int indexToTest =0; indexToTest < asaCalc.getAtomCoords().length; indexToTest++) { + //int indexToTest = 198; + int[] nbsSh = allNbsSh[indexToTest]; + int[] nbs = allNbs[indexToTest]; int countNotInNbs = 0; List listOfMatchingIndices = new ArrayList<>(); for (int i = 0; i < nbsSh.length; i++) { boolean contained = false; for (int j = 0; j < nbs.length; j++) { - if (nbs[j].equals(nbsSh[i])) { + if (nbs[j] == nbsSh[i]) { listOfMatchingIndices.add(j); contained = true; break; @@ -142,7 +144,7 @@ public void testPerformance() throws StructureException, IOException { Structure structure = StructureIO.getStructure("4F5X"); Chain c = structure.getPolyChainByPDB("W"); Atom[] atoms = StructureTools.getAllAtomArray(c); - System.out.printf("Total of %d atoms\n", atoms.length); + System.out.printf("Total of %d atoms. n(n-1)/2= %d \n", atoms.length, atoms.length*(atoms.length-1)/2); int nThreads = 1; // 1. WITH SPATIAL HASHING @@ -164,6 +166,8 @@ public void testPerformance() throws StructureException, IOException { double withSH = totAtoms; System.out.printf("Total ASA is %6.2f \n", totAtoms); + //System.out.println("Distances calculated: " + asaCalc.distancesCalculated); + // 2. WITHOUT SPATIAL HASHING start = System.currentTimeMillis(); @@ -183,6 +187,8 @@ public void testPerformance() throws StructureException, IOException { double withoutSH = totAtoms; System.out.printf("Total ASA is %6.2f \n", totAtoms); + //System.out.println("Distances calculated: " + asaCalc.distancesCalculated); + assertEquals(withoutSH, withSH, 0.000001); } From ddeb7598d4c78893171b0b910e4bd027fa04bbbf Mon Sep 17 00:00:00 2001 From: Jose Duarte Date: Tue, 1 Jan 2019 19:42:00 -0800 Subject: [PATCH 05/11] Efficient neighbor index finding with SH --- .../nbio/structure/asa/AsaCalculator.java | 52 +++++++++++-------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 22b28a47e6..699b44fb82 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -375,34 +375,42 @@ int[][] findNeighborIndices() { */ int[][] findNeighborIndicesSpatialHashing() { - int[][] nbsIndices = new int[atomCoords.length][]; - List contactList = calcContacts(); + Map> indices = new HashMap<>(); + for (Contact contact : contactList) { - for (int k=0; k thisNbIndices = new ArrayList<>(); + List iIndices; + List jIndices; + if (indices.get(i)==null) { + iIndices = new ArrayList<>(); + indices.put(i, iIndices); + } else { + iIndices = indices.get(i); + } + if (indices.get(j)==null) { + jIndices = new ArrayList<>(); + indices.put(j, jIndices); + } else { + jIndices = indices.get(j); + } - // TODO make this the outer loop - for (Contact contact : contactList) { - double dist = contact.getDistance(); - int i; - if (contact.getJ() == k) { - i = contact.getI(); - } else if (contact.getI() == k) { - i = contact.getJ(); - } else { - continue; - } - if (dist < radius + radii[i]) { - thisNbIndices.add(i); - } + double radius = radii[i] + probe + probe; + double dist = contact.getDistance(); + if (dist < radius + radii[j]) { + iIndices.add(j); + jIndices.add(i); } + } - int[] indicesArray = new int[thisNbIndices.size()]; - for (int i=0;i> entry : indices.entrySet()) { + List list = entry.getValue(); + int[] indicesArray = new int[list.size()]; + for (int i=0;i Date: Tue, 1 Jan 2019 23:02:44 -0800 Subject: [PATCH 06/11] Tidying up --- .../biojava/nbio/structure/asa/AsaCalculator.java | 2 +- .../biojava/nbio/structure/asa/TestAsaCalc.java | 15 ++++++--------- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 699b44fb82..8db65b46ae 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -221,7 +221,7 @@ public AsaCalculator(Point3d[] atomCoords, double probe, int nSpherePoints, int */ public GroupAsa[] getGroupAsas() { - TreeMap asas = new TreeMap(); + TreeMap asas = new TreeMap<>(); double[] asasPerAtom = calculateAsas(); diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java index 1f6e4696b4..bc3e1e4842 100644 --- a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java +++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java @@ -141,18 +141,18 @@ public void testPerformance() throws StructureException, IOException { // that's because it depends on the order on how tests were run - JD 2018-03-10 ChemCompGroupFactory.setChemCompProvider(new DownloadChemCompProvider()); - Structure structure = StructureIO.getStructure("4F5X"); - Chain c = structure.getPolyChainByPDB("W"); - Atom[] atoms = StructureTools.getAllAtomArray(c); + Structure structure = StructureIO.getStructure("3HBX"); + Atom[] atoms = StructureTools.getAllAtomArray(structure); System.out.printf("Total of %d atoms. n(n-1)/2= %d \n", atoms.length, atoms.length*(atoms.length-1)/2); int nThreads = 1; - // 1. WITH SPATIAL HASHING + int nSpherePoints = 100; + // 1. WITH SPATIAL HASHING long start = System.currentTimeMillis(); AsaCalculator asaCalc = new AsaCalculator(atoms, AsaCalculator.DEFAULT_PROBE_SIZE, - 100, nThreads); + nSpherePoints, nThreads); asaCalc.setUseSpatialHashingForNeighbors(true); double[] asas = asaCalc.calculateAsas(); @@ -166,14 +166,12 @@ public void testPerformance() throws StructureException, IOException { double withSH = totAtoms; System.out.printf("Total ASA is %6.2f \n", totAtoms); - //System.out.println("Distances calculated: " + asaCalc.distancesCalculated); - // 2. WITHOUT SPATIAL HASHING start = System.currentTimeMillis(); asaCalc = new AsaCalculator(atoms, AsaCalculator.DEFAULT_PROBE_SIZE, - 100, nThreads); + nSpherePoints, nThreads); asaCalc.setUseSpatialHashingForNeighbors(false); asas = asaCalc.calculateAsas(); @@ -187,7 +185,6 @@ public void testPerformance() throws StructureException, IOException { double withoutSH = totAtoms; System.out.printf("Total ASA is %6.2f \n", totAtoms); - //System.out.println("Distances calculated: " + asaCalc.distancesCalculated); assertEquals(withoutSH, withSH, 0.000001); From 041482da4a7b01fb89974d671c8e864f5c85c4ca Mon Sep 17 00:00:00 2001 From: Jose Duarte Date: Wed, 2 Jan 2019 08:44:47 -0800 Subject: [PATCH 07/11] Logging, docs and some minimal optimization --- .../nbio/structure/asa/AsaCalculator.java | 34 ++++++++++++++----- .../nbio/structure/asa/TestAsaCalc.java | 11 +----- 2 files changed, 27 insertions(+), 18 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 8db65b46ae..01c4c9020e 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -251,17 +251,21 @@ public double[] calculateAsas() { double[] asas = new double[atomCoords.length]; if (useSpatialHashingForNeighbors) { + logger.debug("Will use spatial hashing to find neighbors"); neighborIndices = findNeighborIndicesSpatialHashing(); } else { + logger.debug("Will not use spatial hashing to find neighbors"); neighborIndices = findNeighborIndices(); } if (nThreads<=1) { // (i.e. it will also be 1 thread if 0 or negative number specified) + logger.debug("Will use 1 thread for ASA calculation"); for (int i=0;i thisNbIndices = new ArrayList<>(); + List thisNbIndices = new ArrayList<>(initialCapacity); for (int i = 0; i < atomCoords.length; i++) { if (i == k) continue; @@ -372,26 +385,30 @@ int[][] findNeighborIndices() { * Returns the 2-dimensional array with neighbor indices for every atom, * using spatial hashing to avoid all to all distance calculation. * @return 2-dimensional array of size: n_atoms x n_neighbors_per_atom + * @since 5.2.0 */ int[][] findNeighborIndicesSpatialHashing() { + // looking at a typical protein case, number of neighbours are from ~10 to ~50, with an average of ~30 + int initialCapacity = 60; + List contactList = calcContacts(); - Map> indices = new HashMap<>(); + Map> indices = new HashMap<>(atomCoords.length); for (Contact contact : contactList) { - + // note contacts are stored 1-way only, with j>i int i = contact.getI(); int j = contact.getJ(); List iIndices; List jIndices; - if (indices.get(i)==null) { - iIndices = new ArrayList<>(); + if (!indices.containsKey(i)) { + iIndices = new ArrayList<>(initialCapacity); indices.put(i, iIndices); } else { iIndices = indices.get(i); } - if (indices.get(j)==null) { - jIndices = new ArrayList<>(); + if (!indices.containsKey(j)) { + jIndices = new ArrayList<>(initialCapacity); indices.put(j, jIndices); } else { jIndices = indices.get(j); @@ -405,6 +422,7 @@ int[][] findNeighborIndicesSpatialHashing() { } } + // convert map to array for fast access int[][] nbsIndices = new int[atomCoords.length][]; for (Map.Entry> entry : indices.entrySet()) { List list = entry.getValue(); diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java index bc3e1e4842..981820dc08 100644 --- a/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java +++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/asa/TestAsaCalc.java @@ -100,25 +100,16 @@ public void testNeighborIndicesFinding() throws StructureException, IOException int[] nbsSh = allNbsSh[indexToTest]; int[] nbs = allNbs[indexToTest]; - int countNotInNbs = 0; List listOfMatchingIndices = new ArrayList<>(); for (int i = 0; i < nbsSh.length; i++) { - boolean contained = false; for (int j = 0; j < nbs.length; j++) { if (nbs[j] == nbsSh[i]) { listOfMatchingIndices.add(j); - contained = true; break; } } - if (!contained) { - countNotInNbs++; - } } - - //System.out.println("In nbsSh but not in nbs: " + countNotInNbs); - //System.out.println("Number of matching indices: " + listOfMatchingIndices.size()); - + // for (int i = 0; i Date: Wed, 2 Jan 2019 11:11:18 -0800 Subject: [PATCH 08/11] Logging and cleanup --- .../nbio/structure/asa/AsaCalculator.java | 23 ++++++++++++------- .../structure/contact/StructureInterface.java | 4 ++-- .../contact/StructureInterfaceList.java | 20 ++++++++++------ 3 files changed, 30 insertions(+), 17 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index 01c4c9020e..fe7bc230ed 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -134,10 +134,7 @@ public AsaCalculator(Structure structure, double probe, int nSpherePoints, int n radii[i] = getRadius(atoms[i]); } - // initialising the sphere points to sample - spherePoints = generateSpherePoints(nSpherePoints); - - cons = 4.0 * Math.PI / nSpherePoints; + initSpherePoints(nSpherePoints); } /** @@ -169,10 +166,7 @@ public AsaCalculator(Atom[] atoms, double probe, int nSpherePoints, int nThreads radii[i] = getRadius(atoms[i]); } - // initialising the sphere points to sample - spherePoints = generateSpherePoints(nSpherePoints); - - cons = 4.0 * Math.PI / nSpherePoints; + initSpherePoints(nSpherePoints); } /** @@ -206,6 +200,13 @@ public AsaCalculator(Point3d[] atomCoords, double probe, int nSpherePoints, int radii[i] = radius; } + initSpherePoints(nSpherePoints); + } + + private void initSpherePoints(int nSpherePoints) { + + logger.debug("Will use {} sphere points", nSpherePoints); + // initialising the sphere points to sample spherePoints = generateSpherePoints(nSpherePoints); @@ -250,6 +251,7 @@ public double[] calculateAsas() { double[] asas = new double[atomCoords.length]; + long start = System.currentTimeMillis(); if (useSpatialHashingForNeighbors) { logger.debug("Will use spatial hashing to find neighbors"); neighborIndices = findNeighborIndicesSpatialHashing(); @@ -257,7 +259,10 @@ public double[] calculateAsas() { logger.debug("Will not use spatial hashing to find neighbors"); neighborIndices = findNeighborIndices(); } + long end = System.currentTimeMillis(); + logger.debug("Took {} s to find neighbors", (end-start)/1000.0); + start = System.currentTimeMillis(); if (nThreads<=1) { // (i.e. it will also be 1 thread if 0 or negative number specified) logger.debug("Will use 1 thread for ASA calculation"); for (int i=0;i(); - groupAsas2 = new TreeMap(); + groupAsas1 = new TreeMap<>(); + groupAsas2 = new TreeMap<>(); this.totalArea = 0; diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java index 9b4144d96a..b5785b514b 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java @@ -35,7 +35,7 @@ /** * A list of interfaces between 2 molecules (2 sets of atoms) * - * @author duarte_j + * @author Jose Duarte * */ public class StructureInterfaceList implements Serializable, Iterable { @@ -126,8 +126,8 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) { // we get discrepancies (not very big but annoying) which lead to things like negative (small) bsa values - Map uniqAsaChains = new TreeMap(); - Map chainAsas = new TreeMap(); + Map uniqAsaChains = new TreeMap<>(); + Map chainAsas = new TreeMap<>(); // first we gather rotation-unique chains (in terms of AU id and transform id) for (StructureInterface interf:list) { @@ -138,11 +138,15 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) { uniqAsaChains.put(molecId2, interf.getSecondAtomsForAsa(cofactorSizeToUse)); } + logger.debug("Will calculate uncomplexed ASA for {} orientation-unique chains.", uniqAsaChains.size()); + long start = System.currentTimeMillis(); // we only need to calculate ASA for that subset (any translation of those will have same values) for (String molecId:uniqAsaChains.keySet()) { + logger.debug("Calculating uncomplexed ASA for molecId {}, with {} atoms", molecId, uniqAsaChains.get(molecId).length); + AsaCalculator asaCalc = new AsaCalculator(uniqAsaChains.get(molecId), AsaCalculator.DEFAULT_PROBE_SIZE, nSpherePoints, nThreads); @@ -153,8 +157,9 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) { } long end = System.currentTimeMillis(); - logger.debug("Calculated uncomplexed ASA for "+uniqAsaChains.size()+" orientation-unique chains. " - + "Time: "+((end-start)/1000.0)+" s"); + logger.debug("Calculated uncomplexed ASA for {} orientation-unique chains. Time: {} s", uniqAsaChains.size(), ((end-start)/1000.0)); + + logger.debug ("Will calculate complexed ASA for {} pairwise complexes.", list.size()); start = System.currentTimeMillis(); @@ -164,13 +169,14 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) { String molecId1 = interf.getMoleculeIds().getFirst()+interf.getTransforms().getFirst().getTransformId(); String molecId2 = interf.getMoleculeIds().getSecond()+interf.getTransforms().getSecond().getTransformId(); + logger.debug("Calculating complexed ASAs for interface {} between molecules {} and {}", interf.getId(), molecId1, molecId2); + interf.setAsas(chainAsas.get(molecId1), chainAsas.get(molecId2), nSpherePoints, nThreads, cofactorSizeToUse); } end = System.currentTimeMillis(); - logger.debug("Calculated complexes ASA for "+list.size()+" pairwise complexes. " - + "Time: "+((end-start)/1000.0)+" s"); + logger.debug("Calculated complexes ASA for {} pairwise complexes. Time: {} s", list.size(), ((end-start)/1000.0)); // finally we sort based on the ChainInterface.comparable() (based in interfaceArea) From 373a40776c639dde1f01146c68ce92a256bd409a Mon Sep 17 00:00:00 2001 From: Jose Manuel Duarte Date: Wed, 2 Jan 2019 15:04:14 -0800 Subject: [PATCH 09/11] Docs --- .../nbio/structure/asa/AsaCalculator.java | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index fe7bc230ed..c6334a3065 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -54,8 +54,11 @@ public class AsaCalculator { private static final Logger logger = LoggerFactory.getLogger(AsaCalculator.class); - // Bosco uses as default 960, Shrake and Rupley seem to use in their paper 92 (not sure if this is actually the same parameter) - public static final int DEFAULT_N_SPHERE_POINTS = 960; + /** + * The default value for number of sphere points to sample. + * See this paper for a nice study on the effect of this parameter: https://f1000research.com/articles/5-189/v1 + */ + public static final int DEFAULT_N_SPHERE_POINTS = 1000; public static final double DEFAULT_PROBE_SIZE = 1.4; public static final int DEFAULT_NTHREADS = 1; @@ -112,10 +115,11 @@ public void run() { * Constructs a new AsaCalculator. Subsequently call {@link #calculateAsas()} * or {@link #getGroupAsas()} to calculate the ASAs * Only non-Hydrogen atoms are considered in the calculation. - * @param structure - * @param probe - * @param nSpherePoints - * @param nThreads + * @param structure the structure, all non-H atoms will be used + * @param probe the probe size + * @param nSpherePoints the number of points to be used in generating the spherical + * dot-density, the more points the more accurate (and slower) calculation + * @param nThreads the number of parallel threads to use for the calculation * @param hetAtoms if true HET residues are considered, if false they aren't, equivalent to * NACCESS' -h option * @see StructureTools#getAllNonHAtomArray From b674136ec66476cfe5e162d59c29da54ea446405 Mon Sep 17 00:00:00 2001 From: Jose Manuel Duarte Date: Wed, 2 Jan 2019 16:05:50 -0800 Subject: [PATCH 10/11] Docs --- .../org/biojava/nbio/structure/asa/AsaCalculator.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index c6334a3065..b4192775eb 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -37,17 +37,18 @@ /** * Class to calculate Accessible Surface Areas based on * the rolling ball algorithm by Shrake and Rupley. - * + *

* The code is adapted from a python implementation at http://boscoh.com/protein/asapy * (now source is available at https://github.com/boscoh/asa). * Thanks to Bosco K. Ho for a great piece of code and for his fantastic blog. - * + *

* See * Shrake, A., and J. A. Rupley. "Environment and Exposure to Solvent of Protein Atoms. * Lysozyme and Insulin." JMB (1973) 79:351-371. * Lee, B., and Richards, F.M. "The interpretation of Protein Structures: Estimation of * Static Accessibility" JMB (1971) 55:379-400 - * @author duarte_j + * + * @author Jose Duarte * */ public class AsaCalculator { @@ -396,7 +397,6 @@ int[][] findNeighborIndices() { * Returns the 2-dimensional array with neighbor indices for every atom, * using spatial hashing to avoid all to all distance calculation. * @return 2-dimensional array of size: n_atoms x n_neighbors_per_atom - * @since 5.2.0 */ int[][] findNeighborIndicesSpatialHashing() { From 8fb10fdec53cbcb684d0727984b646e602f73918 Mon Sep 17 00:00:00 2001 From: Jose Duarte Date: Thu, 3 Jan 2019 09:07:35 -0800 Subject: [PATCH 11/11] Should be private --- .../main/java/org/biojava/nbio/structure/asa/AsaCalculator.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java index b4192775eb..108cea6067 100644 --- a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java +++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java @@ -63,7 +63,7 @@ public class AsaCalculator { public static final double DEFAULT_PROBE_SIZE = 1.4; public static final int DEFAULT_NTHREADS = 1; - public static final boolean DEFAULT_USE_SPATIAL_HASHING = true; + private static final boolean DEFAULT_USE_SPATIAL_HASHING = true;