Skip to content

Commit 9e501a9

Browse files
josemduartejosemduarte
committed
Further optimization: sort by distance. Achieves ~1/3 of runtime
1 parent c30ab8f commit 9e501a9

File tree

1 file changed

+47
-36
lines changed

1 file changed

+47
-36
lines changed

‎biojava-structure/src/main/java/org/biojava/nbio/structure/asa/AsaCalculator.java‎

Lines changed: 47 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -30,8 +30,7 @@
3030
import java.util.*;
3131
import java.util.concurrent.ExecutorService;
3232
import java.util.concurrent.Executors;
33-
34-
33+
import java.util.stream.Collectors;
3534

3635

3736
/**
@@ -86,8 +85,8 @@ public class AsaCalculator {
8685

8786
private class AsaCalcWorker implements Runnable {
8887

89-
private int i;
90-
private double[] asas;
88+
private final int i;
89+
private final double[] asas;
9190

9291
public AsaCalcWorker(int i, double[] asas) {
9392
this.i = i;
@@ -100,12 +99,21 @@ public void run() {
10099
}
101100
}
102101

102+
private static class IndexAndDistance {
103+
private final int index;
104+
private final double dist;
105+
public IndexAndDistance(int index, double dist) {
106+
this.index = index;
107+
this.dist = dist;
108+
}
109+
}
110+
103111

104-
private Point3d[] atomCoords;
105-
private Atom[] atoms;
106-
private double[] radii;
107-
private double probe;
108-
private int nThreads;
112+
private final Point3d[] atomCoords;
113+
private final Atom[] atoms;
114+
private final double[] radii;
115+
private final double probe;
116+
private final int nThreads;
109117
private Point3d[] spherePoints;
110118
private double cons;
111119
private int[][] neighborIndices;
@@ -123,7 +131,6 @@ public void run() {
123131
* @param nThreads the number of parallel threads to use for the calculation
124132
* @param hetAtoms if true HET residues are considered, if false they aren't, equivalent to
125133
* NACCESS' -h option
126-
* @see StructureTools#getAllNonHAtomArray
127134
*/
128135
public AsaCalculator(Structure structure, double probe, int nSpherePoints, int nThreads, boolean hetAtoms) {
129136
this.atoms = StructureTools.getAllNonHAtomArray(structure, hetAtoms);
@@ -312,10 +319,8 @@ public double[] calculateAsas() {
312319
//12 threads, time: 0.9s -- x11.4
313320

314321

315-
316322
ExecutorService threadPool = Executors.newFixedThreadPool(nThreads);
317323

318-
319324
for (int i=0;i<atomCoords.length;i++) {
320325
threadPool.submit(new AsaCalcWorker(i,asas));
321326
}
@@ -404,14 +409,14 @@ int[][] findNeighborIndicesSpatialHashing() {
404409
int initialCapacity = 60;
405410

406411
List<Contact> contactList = calcContacts();
407-
Map<Integer, List<Integer>> indices = new HashMap<>(atomCoords.length);
412+
Map<Integer, List<IndexAndDistance>> indices = new HashMap<>(atomCoords.length);
408413
for (Contact contact : contactList) {
409414
// note contacts are stored 1-way only, with j>i
410415
int i = contact.getI();
411416
int j = contact.getJ();
412417

413-
List<Integer> iIndices;
414-
List<Integer> jIndices;
418+
List<IndexAndDistance> iIndices;
419+
List<IndexAndDistance> jIndices;
415420
if (!indices.containsKey(i)) {
416421
iIndices = new ArrayList<>(initialCapacity);
417422
indices.put(i, iIndices);
@@ -428,17 +433,21 @@ int[][] findNeighborIndicesSpatialHashing() {
428433
double radius = radii[i] + probe + probe;
429434
double dist = contact.getDistance();
430435
if (dist < radius + radii[j]) {
431-
iIndices.add(j);
432-
jIndices.add(i);
436+
iIndices.add(new IndexAndDistance(j, dist));
437+
jIndices.add(new IndexAndDistance(i, dist));
433438
}
434439
}
435440

436441
// convert map to array for fast access
437442
int[][] nbsIndices = new int[atomCoords.length][];
438-
for (Map.Entry<Integer, List<Integer>> entry : indices.entrySet()) {
439-
List<Integer> list = entry.getValue();
443+
for (Map.Entry<Integer, List<IndexAndDistance>> entry : indices.entrySet()) {
444+
List<IndexAndDistance> list = entry.getValue();
445+
// Sorting by closest to farthest away neighbors achieves faster runtimes when checking for occluded sphere sample points in calcSingleAsa.
446+
// This follows the ideas exposed in Eisenhaber et al, J Comp Chemistry 1994 (https://onlinelibrary.wiley.com/doi/epdf/10.1002/jcc.540160303)
447+
list = list.stream().sorted(Comparator.comparingDouble(o -> o.dist)).collect(Collectors.toList());
440448
int[] indicesArray = new int[list.size()];
441-
for (int i=0;i<entry.getValue().size();i++) indicesArray[i] = list.get(i);
449+
for (int i=0; i<list.size(); i++) indicesArray[i] = list.get(i).index;
450+
442451
nbsIndices[entry.getKey()] = indicesArray;
443452
}
444453

@@ -475,34 +484,29 @@ private double calcSingleAsa(int i) {
475484

476485
int n_neighbor = neighborIndices[i].length;
477486
int[] neighbor_indices = neighborIndices[i];
478-
int j_closest_neighbor = 0;
479487
double radius = probe + radii[i];
480488

481489
int n_accessible_point = 0;
482490

491+
int[] numDistsCalced = null;
492+
if (logger.isDebugEnabled()) numDistsCalced = new int[n_neighbor];
493+
483494
for (Point3d point: spherePoints){
484495
boolean is_accessible = true;
485496
Point3d test_point = new Point3d(point.x*radius + atom_i.x,
486497
point.y*radius + atom_i.y,
487498
point.z*radius + atom_i.z);
488499

489-
int[] cycled_indices = new int[n_neighbor];
490-
int arind = 0;
491-
for (int ind=j_closest_neighbor;ind<n_neighbor;ind++) {
492-
cycled_indices[arind] = ind;
493-
arind++;
494-
}
495-
for (int ind=0;ind<j_closest_neighbor;ind++){
496-
cycled_indices[arind] = ind;
497-
arind++;
498-
}
499-
500-
for (int j: cycled_indices) {
501-
Point3d atom_j = atomCoords[neighbor_indices[j]];
502-
double r = radii[neighbor_indices[j]] + probe;
500+
int iter = -1;
501+
for (int j : neighbor_indices) {
502+
iter++;
503+
Point3d atom_j = atomCoords[j];
504+
double r = radii[j] + probe;
503505
double diff_sq = test_point.distanceSquared(atom_j);
506+
507+
if (numDistsCalced!=null) numDistsCalced[iter]++;
508+
504509
if (diff_sq < r*r) {
505-
j_closest_neighbor = j;
506510
is_accessible = false;
507511
break;
508512
}
@@ -511,6 +515,13 @@ private double calcSingleAsa(int i) {
511515
n_accessible_point++;
512516
}
513517
}
518+
519+
if (numDistsCalced!=null) {
520+
int sum = 0;
521+
for (int j = 0; j < n_neighbor; j++) sum += numDistsCalced[j];
522+
logger.debug("Number of sample points distances calculated for neighbors of i={} : average {}, all {}", i, (double) sum / (double) n_neighbor, numDistsCalced);
523+
}
524+
514525
return cons*n_accessible_point*radius*radius;
515526
}
516527

0 commit comments

Comments
 (0)