Merge pull request #964 from aidenlab/mvm2

Mvm2

src/juicebox/HiC.java

@@ -1235,10 +1235,7 @@ public void clearAllMatrixZoomDataCache() {
     }
 
     private void clearAllCacheForDataset(Dataset ds) {
-        Matrix matrix = ds.getMatrix(xContext.getChromosome(), yContext.getChromosome());
-        for (HiCZoom zoom : ds.getBpZooms()) {
-            matrix.getZoomData(zoom).clearCache();
-        }
+        ds.clearCache(false);
     }
 
     public List<Pair<MotifAnchor, MotifAnchor>> getRTreeHandlerIntersectingFeatures(String name, int g1, int g2) {

src/juicebox/HiCGlobals.java

@@ -36,7 +36,7 @@
  */
 public class HiCGlobals {
 
-    public static final String versionNum = "2.11.00";
+    public static final String versionNum = "2.12.00";
     public static final String juiceboxTitle = "[Juicebox " + versionNum + "] Hi-C Map ";
 
     // MainWindow variables

src/juicebox/data/Dataset.java

@@ -1003,4 +1003,26 @@ public NormalizationHandler getNormalizationHandler() {
     public int getDepthBase() {
         return v9DepthBase;
     }
+
+    public void clearCache(boolean onlyClearInter) {
+        for (Matrix matrix : matrices.values()) {
+            for (HiCZoom zoom : getBpZooms()) {
+                try {
+                    matrix.getZoomData(zoom).clearCache(onlyClearInter);
+                } catch (Exception e) {
+                    System.err.println("Clearing err: " + e.getLocalizedMessage());
+                }
+            }
+        }
+    }
+
+    public void clearCache(boolean onlyClearInter, HiCZoom zoom) {
+        for (Matrix matrix : matrices.values()) {
+            try {
+                matrix.getZoomData(zoom).clearCache(onlyClearInter);
+            } catch (Exception e) {
+                System.err.println("Clearing z_err: " + e.getLocalizedMessage());
+            }
+        }
+    }
 }

src/juicebox/data/DatasetReaderV2.java

@@ -63,12 +63,12 @@ public class DatasetReaderV2 extends AbstractDatasetReader {
      * Cache of chromosome name -> array of restriction sites
      */
     private final Map<String, int[]> fragmentSitesCache = new HashMap<>();
-    private final Map<String, IndexEntry> masterIndex = Collections.synchronizedMap(new HashMap<>());
+    private final Map<String, IndexEntry> masterIndex = new HashMap<>();
     private Map<String, LargeIndexEntry> normVectorIndex;
     private final Dataset dataset;
     private int version = -1;
     private Map<String, FragIndexEntry> fragmentSitesIndex;
-    private final Map<String, BlockIndex> blockIndexMap = Collections.synchronizedMap(new HashMap<>());
+    private final Map<String, BlockIndex> blockIndexMap = new HashMap<>();
     private long masterIndexPos;
     private long normVectorFilePosition;
     private long nviHeaderPosition;
@@ -332,11 +332,15 @@ private Pair<MatrixZoomData, Long> readMatrixZoomData(Chromosome chr1, Chromosom
         if (binSize < 50 && HiCGlobals.allowDynamicBlockIndex) {
             int maxPossibleBlockNumber = blockColumnCount * blockColumnCount - 1;
             DynamicBlockIndex blockIndex = new DynamicBlockIndex(getValidStream(), nBlocks, maxPossibleBlockNumber, currentFilePointer);
-            blockIndexMap.put(zd.getKey(), blockIndex);
+            synchronized (blockIndexMap) {
+                blockIndexMap.put(zd.getKey(), blockIndex);
+            }
         } else {
             BlockIndex blockIndex = new BlockIndex(nBlocks);
             blockIndex.populateBlocks(dis);
-            blockIndexMap.put(zd.getKey(), blockIndex);
+            synchronized (blockIndexMap) {
+                blockIndexMap.put(zd.getKey(), blockIndex);
+            }
         }
         currentFilePointer += (nBlocks * 16L);
 
@@ -483,7 +487,9 @@ private void readFooter(long position) throws IOException {
             long filePosition = dis.readLong();
             int sizeInBytes = dis.readInt();
             currentPosition += 12;
-            masterIndex.put(key, new IndexEntry(filePosition, sizeInBytes));
+            synchronized (masterIndex) {
+                masterIndex.put(key, new IndexEntry(filePosition, sizeInBytes));
+            }
         }
 
         Map<String, ExpectedValueFunction> expectedValuesMap = new LinkedHashMap<>();

src/juicebox/data/MatrixZoomData.java

@@ -34,6 +34,7 @@
 import juicebox.data.basics.Chromosome;
 import juicebox.data.iterator.IteratorContainer;
 import juicebox.data.iterator.ListOfListGenerator;
+import juicebox.data.iterator.ZDIteratorContainer;
 import juicebox.data.v9depth.LogDepth;
 import juicebox.data.v9depth.V9Depth;
 import juicebox.gui.SuperAdapter;
@@ -61,11 +62,6 @@
 import java.util.concurrent.atomic.AtomicInteger;
 
 
-
-/**
- * @author jrobinso
- * @since Aug 10, 2010
- */
 public class MatrixZoomData {
 
     final Chromosome chr1;  // Chromosome on the X axis
@@ -1197,8 +1193,15 @@ public void setAverageCount(double averageCount) {
         this.averageCount = averageCount;
     }
 
-    public void clearCache() {
-        blockCache.clear();
+    public void clearCache(boolean onlyClearInter) {
+        if (onlyClearInter && isIntra) return;
+        if (HiCGlobals.useCache) {
+            blockCache.clear();
+        }
+        if (iteratorContainer != null) {
+            iteratorContainer.clear();
+            iteratorContainer = null;
+        }
     }
 
     private Iterator<ContactRecord> getNewContactRecordIterator() {
@@ -1212,4 +1215,8 @@ public IteratorContainer getIteratorContainer() {
         }
         return iteratorContainer;
     }
+
+    public IteratorContainer getFromFileIteratorContainer() {
+        return new ZDIteratorContainer(reader, this, blockCache);
+    }
 }

src/juicebox/data/basics/ListOfFloatArrays.java

@@ -151,8 +151,10 @@ public void addTo(long index, float value) {
     public void addValuesFrom(ListOfFloatArrays other) {
         if (overallLength == other.overallLength) {
             for (int i = 0; i < internalList.size(); i++) {
-                for (int j = 0; j < internalList.get(i).length; j++) {
-                    internalList.get(i)[j] += other.internalList.get(i)[j];
+                float[] array = internalList.get(i);
+                float[] otherArray = other.internalList.get(i);
+                for (int j = 0; j < array.length; j++) {
+                    array[j] += otherArray[j];
                 }
             }
         } else {
@@ -184,8 +186,12 @@ public void multiplyEverythingBy(double val) {
     public ListOfDoubleArrays convertToDoubles() {
         ListOfDoubleArrays newList = new ListOfDoubleArrays(overallLength);
         for (int j = 0; j < internalList.size(); j++) {
-            for (int k = 0; k < internalList.get(j).length; k++) {
-                newList.getValues().get(j)[k] = internalList.get(j)[k];
+
+            float[] array = internalList.get(j);
+            double[] newArray = newList.getValues().get(j);
+
+            for (int k = 0; k < array.length; k++) {
+                newArray[k] = array[k];
             }
         }
         return newList;

src/juicebox/data/iterator/CoupledIteratorAndOffset.java

@@ -0,0 +1,54 @@
+/*
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2011-2021 Broad Institute, Aiden Lab, Rice University, Baylor College of Medicine
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ *  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ *  FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT. IN NO EVENT SHALL THE
+ *  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ *  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ *  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ *  THE SOFTWARE.
+ */
+
+package juicebox.data.iterator;
+
+import juicebox.data.ContactRecord;
+
+import java.util.Iterator;
+
+public class CoupledIteratorAndOffset implements Iterator<ContactRecord> {
+
+    private final Iterator<ContactRecord> internalIterator;
+    private final int xOffset, yOffset;
+
+    public CoupledIteratorAndOffset(Iterator<ContactRecord> iterator, int xOffset, int yOffset) {
+        internalIterator = iterator;
+        this.xOffset = xOffset;
+        this.yOffset = yOffset;
+    }
+
+    @Override
+    public boolean hasNext() {
+        return internalIterator.hasNext();
+    }
+
+    @Override
+    public ContactRecord next() {
+        ContactRecord cr = internalIterator.next();
+        int binX = cr.getBinX() + xOffset;
+        int binY = cr.getBinY() + yOffset;
+        return new ContactRecord(binX, binY, cr.getCounts());
+    }
+}

src/juicebox/data/iterator/GWIteratorContainer.java

@@ -28,9 +28,13 @@
 import juicebox.data.ContactRecord;
 import juicebox.data.Dataset;
 import juicebox.data.basics.Chromosome;
+import juicebox.data.basics.ListOfFloatArrays;
+import juicebox.tools.dev.ParallelizedJuicerTools;
 import juicebox.windowui.HiCZoom;
 
 import java.util.Iterator;
+import java.util.List;
+import java.util.concurrent.atomic.AtomicInteger;
 
 
 public class GWIteratorContainer extends IteratorContainer {
@@ -62,4 +66,38 @@ private static long calculateMatrixSize(ChromosomeHandler handler, HiCZoom zoom)
     public Iterator<ContactRecord> getNewContactRecordIterator() {
         return new GenomeWideIterator(dataset, handler, zoom, includeIntra);
     }
+
+    public List<Iterator<ContactRecord>> getAllFromFileContactRecordIterators() {
+        return GenomeWideIterator.getAllFromFileIterators(dataset, handler, zoom, includeIntra);
+    }
+
+    @Override
+    public ListOfFloatArrays sparseMultiply(ListOfFloatArrays vector, long vectorLength) {
+        final ListOfFloatArrays totalSumVector = new ListOfFloatArrays(vectorLength);
+
+        List<Iterator<ContactRecord>> allIterators = getAllFromFileContactRecordIterators();
+
+        AtomicInteger index = new AtomicInteger(0);
+        ParallelizedJuicerTools.launchParallelizedCode(numCPUMatrixThreads, () -> {
+            int i = index.getAndIncrement();
+            ListOfFloatArrays accumSumVector = new ListOfFloatArrays(vectorLength);
+            while (i < allIterators.size()) {
+                accumSumVector.addValuesFrom(ZDIteratorContainer.matrixVectorMultiplyOnIterator(
+                        allIterators.get(i), vector, vectorLength));
+                i = index.getAndIncrement();
+            }
+            synchronized (totalSumVector) {
+                totalSumVector.addValuesFrom(accumSumVector);
+            }
+        });
+
+        allIterators.clear();
+
+        return totalSumVector;
+    }
+
+    @Override
+    public void clear() {
+        // null, doesn't need to clean anything
+    }
 }

src/juicebox/data/iterator/GenomeWideIterator.java

@@ -28,7 +28,9 @@
 import juicebox.data.basics.Chromosome;
 import juicebox.windowui.HiCZoom;
 
+import java.util.ArrayList;
 import java.util.Iterator;
+import java.util.List;
 
 public class GenomeWideIterator implements Iterator<ContactRecord> {
 
@@ -62,6 +64,35 @@ public boolean hasNext() {
         return getNextIterator();
     }
 
+    public static List<Iterator<ContactRecord>> getAllFromFileIterators(Dataset dataset, ChromosomeHandler handler,
+                                                                        HiCZoom zoom, boolean includeIntra) {
+        Chromosome[] chromosomes = handler.getChromosomeArrayWithoutAllByAll();
+        List<Iterator<ContactRecord>> allIterators = new ArrayList<>();
+
+        int xOffset = 0;
+        for (int i = 0; i < chromosomes.length; i++) {
+            Chromosome c1 = chromosomes[i];
+            int yOffset = 0 + xOffset;
+            for (int j = i; j < chromosomes.length; j++) {
+                Chromosome c2 = chromosomes[j];
+
+                if (c1.getIndex() < c2.getIndex() || (c1.equals(c2) && includeIntra)) {
+                    MatrixZoomData zd = HiCFileTools.getMatrixZoomData(dataset, c1, c2, zoom);
+                    if (zd != null) {
+                        IteratorContainer ic = zd.getFromFileIteratorContainer();
+                        Iterator<ContactRecord> iterator = ic.getNewContactRecordIterator();
+                        if (iterator != null && iterator.hasNext()) {
+                            allIterators.add(new CoupledIteratorAndOffset(iterator, xOffset, yOffset));
+                        }
+                    }
+                }
+                yOffset += c2.getLength() / zoom.getBinSize() + 1;
+            }
+            xOffset += c1.getLength() / zoom.getBinSize() + 1;
+        }
+        return allIterators;
+    }
+
     private boolean getNextIterator() {
         while (c1i < chromosomes.length) {
             Chromosome c1 = chromosomes[c1i];
@@ -71,8 +102,9 @@ private boolean getNextIterator() {
                 if (c1.getIndex() < c2.getIndex() || (c1.equals(c2) && includeIntra)) {
                     MatrixZoomData zd = HiCFileTools.getMatrixZoomData(dataset, c1, c2, zoom);
                     if (zd != null) {
-                        currentIterator = zd.getIteratorContainer().getNewContactRecordIterator();
-                        if (currentIterator != null && currentIterator.hasNext()) {
+                        Iterator<ContactRecord> newIterator = zd.getFromFileIteratorContainer().getNewContactRecordIterator();
+                        if (newIterator != null && newIterator.hasNext()) {
+                            currentIterator = new CoupledIteratorAndOffset(newIterator, recentAddX, recentAddY);
                             return true;
                         }
                     }
@@ -81,18 +113,15 @@ private boolean getNextIterator() {
                 c2i++;
             }
             recentAddX += c1.getLength() / zoom.getBinSize() + 1;
-            recentAddY = 0;
+            recentAddY = 0 + recentAddX;
             c1i++;
-            c2i = 0;
+            c2i = c1i;
         }
         return false;
     }
 
     @Override
     public ContactRecord next() {
-        ContactRecord cr = currentIterator.next();
-        int binX = cr.getBinX() + recentAddX;
-        int binY = cr.getBinY() + recentAddY;
-        return new ContactRecord(binX, binY, cr.getCounts());
+        return currentIterator.next();
     }
 }

src/juicebox/data/iterator/IteratorContainer.java

@@ -25,13 +25,15 @@
 package juicebox.data.iterator;
 
 import juicebox.data.ContactRecord;
+import juicebox.data.basics.ListOfFloatArrays;
 
 import java.util.Iterator;
 
 public abstract class IteratorContainer {
 
     private final long matrixSize;
     private long numberOfContactRecords = -1;
+    public static int numCPUMatrixThreads = 10;
 
     public IteratorContainer(long matrixSize) {
         this.matrixSize = matrixSize;
@@ -66,4 +68,9 @@ public boolean getIsThereEnoughMemoryForNormCalculation() {
         // float is 4 bytes; one for each row
         return matrixSize * 4 < Runtime.getRuntime().maxMemory();
     }
+
+    public abstract ListOfFloatArrays sparseMultiply(ListOfFloatArrays vector, long vectorLength);
+
+    public abstract void clear();
+
 }