The label is used when generating a diagram to illustrate the search space and costs
+ */
+ public Set> subFinderLabels();
+
/**
* Returns the total number of possible sub-haplotypes.
* @return 0 or greater.
@@ -67,5 +94,22 @@ interface KBestSubHaplotypeFinder {
*
* @return never {@code null}.
*/
- public abstract KBestHaplotype getKBest(int k);
+ public KBestHaplotype getKBest(int k);
+
+ /**
+ * Checks whether the top vertex for this finder is a reference haplotype vertex.
+ *
+ * @return {@code true} iff the top vertex for this finder is a reference vertex.
+ */
+ public boolean isReference();
+
+ /**
+ * Calculate the score of a sequence of bases.
+ *
+ * @param bases array containing the query base sequence.
+ * @param offset first position of the query base sequence in {@code bases} .
+ * @param length length of the query base sequence.
+ * @return {@link Double#NaN} if there is no score for this sequence, otherwise a valid score value.
+ */
+ public double score(byte[] bases, int offset, int length);
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/MultiSampleEdge.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/MultiSampleEdge.java
index 978d83eb4..657ecfd85 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/MultiSampleEdge.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/MultiSampleEdge.java
@@ -49,20 +49,24 @@ package org.broadinstitute.sting.gatk.walkers.haplotypecaller.graphs;
import java.util.PriorityQueue;
/**
- * edge class for connecting nodes in the graph that tracks some per-sample information
- *
+ * Edge class for connecting nodes in the graph that tracks some per-sample information.
+ *
* This class extends BaseEdge with the additional functionality of tracking the maximum
* multiplicity seen within any single sample. The workflow for using this class is:
- *
- * MultiSampleEdge e = new MultiSampleEdge(ref, 1)
- * e.incMultiplicity(1) // total is 2, per sample is 2, max per sample is 1
- * e.getPruningMultiplicity() // = 1
- * e.flushSingleSampleMultiplicity() // total is 2, per sample is 0, max per sample is 2
- * e.getPruningMultiplicity() // = 2
- * e.incMultiplicity(3) // total is 5, per sample is 3, max per sample is 2
- * e.getPruningMultiplicity() // = 2
- * e.flushSingleSampleMultiplicity() // total is 5, per sample is 0, max per sample is 3
- * e.getPruningMultiplicity() // = 3
+ *
+ *
+ * {@code
+ * MultiSampleEdge e = new MultiSampleEdge(ref, 1)
+ * e.incMultiplicity(1) // total is 2, per sample is 2, max per sample is 1
+ * e.getPruningMultiplicity() // = 1
+ * e.flushSingleSampleMultiplicity() // total is 2, per sample is 0, max per sample is 2
+ * e.getPruningMultiplicity() // = 2
+ * e.incMultiplicity(3) // total is 5, per sample is 3, max per sample is 2
+ * e.getPruningMultiplicity() // = 2
+ * e.flushSingleSampleMultiplicity() // total is 5, per sample is 0, max per sample is 3
+ * e.getPruningMultiplicity() // = 3
+ * }
+ *
*/
public class MultiSampleEdge extends BaseEdge {
private int currentSingleSampleMultiplicity;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/RecursiveSubHaplotypeFinder.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/RecursiveSubHaplotypeFinder.java
index 0fbbfdc64..fe85befef 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/RecursiveSubHaplotypeFinder.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/RecursiveSubHaplotypeFinder.java
@@ -45,9 +45,10 @@
*/
package org.broadinstitute.sting.gatk.walkers.haplotypecaller.graphs;
-import java.util.ArrayList;
-import java.util.Collection;
-import java.util.Map;
+import org.broadinstitute.sting.utils.Utils;
+import org.broadinstitute.sting.utils.collections.Pair;
+
+import java.util.*;
/**
* General recursive sub-haplotype finder.
@@ -67,7 +68,11 @@ import java.util.Map;
*
* @author Valentin Ruano-Rubio <valentin@broadinstitute.org>
*/
-class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
+class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
+
+
+ private final SeqVertex vertex;
+ private final boolean isReference;
/**
* Creates a recursive sub-haplotype finder give the target graph, first vertex and all possible outgoing edges
@@ -80,20 +85,83 @@ class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
* @param vertex first vertex for all sub-haplotype solutions provided by this finder
* @param children map from outgoing edge to the corresponding sub-sub-haplotype finder.
*/
- public RecursiveSubHaplotypeFinder(final SeqVertex vertex,
+ public RecursiveSubHaplotypeFinder(final SeqGraph graph, final SeqVertex vertex,
final Map children) {
super(createChildFinderCollection(vertex, children));
+ this.vertex = vertex;
+ this.isReference = graph.isReferenceNode(vertex);
}
- private static Collection createChildFinderCollection(final SeqVertex vertex, final Map finders) {
+ /**
+ * Wraps the descendant vertices finders in order to take advantage of the {@link AggregatedSubHaplotypeFinder}
+ * common code.
+ *
+ * Automatically calibrates the edge score (cost) so that it takes into account the total across all edges.
+ *
+ *
+ * @param vertex the parent vertex.
+ * @param finders the child vertices indexed by the connecting edge.
+ * @return never {@code null} but potentially an empty collection if there is child returning some sub-haplotype
+ * solution.
+ */
+ private static Collection createChildFinderCollection(final SeqVertex vertex,
+ final Map finders) {
if (finders == null) throw new IllegalArgumentException("the edge to child map cannot be null");
- final Collection result = new ArrayList<>(finders.size());
- for (final Map.Entry e : finders.entrySet())
- result.add(new EdgeSubHaplotypeFinder(vertex,e.getKey(), e.getValue()));
+ final ArrayList result = new ArrayList<>(finders.size());
+ for (final Map.Entry e : finders.entrySet()) {
+ final EdgeSubHaplotypeFinder subFinder = new EdgeSubHaplotypeFinder(vertex,e.getKey(), e.getValue());
+ if (subFinder.getCount() == 0) continue;
+ result.add(subFinder);
+ }
+ if (result.size() == 0)
+ return Collections.emptySet();
+ else if (result.size() == 1) // no calibration needed, by default edgeScore is 0.
+ return Collections.singleton(result.get(0));
+ else {
+ double totalEdgeMultiplicityAcrossEdges = 0;
+ for (final EdgeSubHaplotypeFinder finder : result)
+ totalEdgeMultiplicityAcrossEdges += Math.max(0.5,finder.edge.getMultiplicity());
+ final double log10TotalEdgeMultiplicityAcrossEdges = Math.log10(totalEdgeMultiplicityAcrossEdges);
+ for (final EdgeSubHaplotypeFinder finder : result)
+ finder.calibrateEdgeScore(log10TotalEdgeMultiplicityAcrossEdges);
+ return result;
+ }
+ }
+
+ @Override
+ public boolean isReference() {
+ return isReference;
+ }
+
+ @Override
+ public String label() {
+ return vertex.getSequenceString();
+ }
+
+ @Override
+ public Set> subFinderLabels() {
+ final Set> result = new LinkedHashSet<>(subFinders.size());
+ for (final EdgeSubHaplotypeFinder subFinder : subFinders)
+ result.add(new Pair<>(subFinder,simplifyZeros(String.format("%.4f", subFinder.edgeScore))));
return result;
}
- private static class EdgeSubHaplotypeFinder implements KBestSubHaplotypeFinder {
+ /**
+ * Removes zeros decimal positions from edge-labels.
+ *
+ * @param edgeLabel the original label to reformat.
+ * @return never {@code null}, the reformatted label.
+ */
+ private String simplifyZeros(final String edgeLabel) {
+ if (edgeLabel.equals("0.000") || edgeLabel.equals("-0.000") )
+ return "0.";
+ int i = edgeLabel.length() - 1;
+ while (edgeLabel.charAt(i) == '0')
+ i--;
+ return (i == edgeLabel.length() - 1) ? edgeLabel : edgeLabel.substring(0,i);
+ }
+
+ protected static class EdgeSubHaplotypeFinder implements KBestSubHaplotypeFinder {
private final KBestSubHaplotypeFinder childFinder;
@@ -101,10 +169,32 @@ class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
private final BaseEdge edge;
+ private double edgeScore = 0;
+
private EdgeSubHaplotypeFinder(final SeqVertex vertex, final BaseEdge edge, final KBestSubHaplotypeFinder childFinder) {
this.childFinder = childFinder;
this.edge = edge;
this.vertex = vertex;
+ this.edgeScore = 0;
+ }
+
+ private void calibrateEdgeScore(final double log10TotalMultiplicityAcrossOutgoingEdges) {
+ edgeScore = Math.log10(Math.max(edge.getMultiplicity(),0.5)) - log10TotalMultiplicityAcrossOutgoingEdges;
+ }
+
+ @Override
+ public String id() {
+ return childFinder.id();
+ }
+
+ @Override
+ public String label() {
+ return childFinder.label();
+ }
+
+ @Override
+ public Set> subFinderLabels() {
+ return childFinder.subFinderLabels();
}
@Override
@@ -114,8 +204,31 @@ class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
@Override
public KBestHaplotype getKBest(int k) {
- return new ChildKBestSubHaplotype(vertex,edge,childFinder.getKBest(k));
+ return new ChildKBestSubHaplotype(vertex,edge,childFinder.getKBest(k),edgeScore);
}
+
+ @Override
+ public boolean isReference() {
+ return childFinder.isReference();
+ }
+
+ @Override
+ public double score(final byte[] bases, final int offset, final int length) {
+ if (length == 0)
+ return 0;
+ final byte[] vertexSequence = vertex.getSequence();
+ if (length < vertexSequence.length) // query is not long enough to have any score.
+ return Double.NaN;
+ else if (!Utils.equalRange(vertexSequence,0,bases,offset,vertexSequence.length))
+ return Double.NaN;
+ else
+ return edgeScore + childFinder.score(bases,offset + vertexSequence.length,length - vertexSequence.length);
+ }
+ }
+
+ @Override
+ public String id() {
+ return "v" + vertex.getId();
}
/**
@@ -129,13 +242,14 @@ class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
*/
private static class ChildKBestSubHaplotype extends KBestHaplotype {
- private final int score;
+ private final double score;
private final KBestHaplotype child;
private final SeqVertex vertex;
private final boolean isReference;
- public ChildKBestSubHaplotype(final SeqVertex vertex, final BaseEdge edge, final KBestHaplotype child) {
- this.score = edge.getMultiplicity() + child.score();
+
+ public ChildKBestSubHaplotype(final SeqVertex vertex, final BaseEdge edge, final KBestHaplotype child, final double edgeScore) {
+ this.score = edgeScore + child.score();
this.vertex = vertex;
this.child = child;
this.isReference = edge.isRef() && child.isReference();
@@ -147,7 +261,7 @@ class RecursiveSubHaplotypeFinder extends AggregatedSubHaplotypeFinder {
}
@Override
- public int score() {
+ public double score() {
return score;
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SeqGraph.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SeqGraph.java
index c8c6abb86..087c07be7 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SeqGraph.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SeqGraph.java
@@ -291,16 +291,9 @@ public class SeqGraph extends BaseGraph {
final SeqVertex addedVertex = mergeLinearChainVertices(linearChain);
addVertex(addedVertex);
- final Set inEdges = incomingEdgesOf(first);
- final Set outEdges = outgoingEdgesOf(last);
-
- final int nEdges = inEdges.size() + outEdges.size();
- int sharedWeightAmongEdges = nEdges == 0 ? 0 : sumEdgeWeightAlongChain(linearChain) / nEdges;
- final BaseEdge inc = new BaseEdge(false, sharedWeightAmongEdges); // template to make .add function call easy
-
// update the incoming and outgoing edges to point to the new vertex
- for( final BaseEdge edge : outEdges ) { addEdge(addedVertex, getEdgeTarget(edge), edge.copy().add(inc)); }
- for( final BaseEdge edge : inEdges ) { addEdge(getEdgeSource(edge), addedVertex, edge.copy().add(inc)); }
+ for( final BaseEdge edge : outgoingEdgesOf(last) ) { addEdge(addedVertex, getEdgeTarget(edge), edge.copy()); }
+ for( final BaseEdge edge : incomingEdgesOf(first) ) { addEdge(getEdgeSource(edge), addedVertex, edge.copy()); }
removeAllVertices(linearChain);
return true;
@@ -313,29 +306,6 @@ public class SeqGraph extends BaseGraph {
return new SeqVertex( seqsCat );
}
- /**
- * Get the sum of the edge weights on a linear chain of at least 2 elements
- *
- * @param chain a linear chain of vertices with at least 2 vertices
- * @return the sum of the multiplicities along all edges connecting vertices within the chain
- */
- @Requires({"chain != null", "chain.size() >= 2"})
- private int sumEdgeWeightAlongChain(final LinkedList chain) {
- int sum = 0;
- SeqVertex prev = null;
-
- for ( final SeqVertex v : chain ) {
- if ( prev != null ) {
- final BaseEdge e = getEdge(prev, v);
- if ( e == null ) throw new IllegalStateException("Something wrong with the linear chain, got a null edge between " + prev + " and " + v);
- sum += e.getMultiplicity();
- }
- prev = v;
- }
-
- return sum;
- }
-
/**
* Base class for transformation operations that need to iterate over proposed vertices, where
* each proposed vertex is a seed vertex for a potential transformation.
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitter.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitter.java
index 284062749..401cbf18c 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitter.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitter.java
@@ -247,12 +247,12 @@ public class SharedVertexSequenceSplitter {
if ( needPrefixNode ) {
outer.addVertex(prefixV);
- if ( top != null ) outer.addEdge(top, prefixV, BaseEdge.orRef(splitGraph.outgoingEdgesOf(prefixV), 0));
+ if ( top != null ) outer.addEdge(top, prefixV, BaseEdge.orRef(splitGraph.outgoingEdgesOf(prefixV), 1));
}
if ( needSuffixNode ) {
outer.addVertex(suffixV);
- if ( bot != null ) outer.addEdge(suffixV, bot, BaseEdge.orRef(splitGraph.incomingEdgesOf(suffixV), 0));
+ if ( bot != null ) outer.addEdge(suffixV, bot, BaseEdge.orRef(splitGraph.incomingEdgesOf(suffixV), 1));
}
if ( topForConnect != null ) {
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java
index a7989ac2c..e03e26e0a 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java
@@ -52,7 +52,6 @@ import org.broadinstitute.sting.gatk.walkers.haplotypecaller.Kmer;
import org.broadinstitute.sting.gatk.walkers.haplotypecaller.graphs.*;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.jgrapht.alg.CycleDetector;
import java.io.File;
import java.util.*;
@@ -88,8 +87,7 @@ public class ReadThreadingGraph extends DanglingChainMergingGraph implements Kme
/**
*
*/
-
- final boolean debugGraphTransformations;
+ private final boolean debugGraphTransformations;
final byte minBaseQualityToUseInAssembly;
protected boolean increaseCountsBackwards = true;
@@ -319,13 +317,6 @@ public class ReadThreadingGraph extends DanglingChainMergingGraph implements Kme
removeAllVertices(verticesToRemove);
}
- /**
- * @return true if the graph has cycles, false otherwise
- */
- public boolean hasCycles() {
- return new CycleDetector<>(this).detectCycles();
- }
-
/**
* Does the graph not have enough complexity? We define low complexity as a situation where the number
* of non-unique kmers is more than 20% of the total number of kmers.
@@ -419,39 +410,10 @@ public class ReadThreadingGraph extends DanglingChainMergingGraph implements Kme
return counter.getKmersWithCountsAtLeast(2);
}
- /**
- * Convert this kmer graph to a simple sequence graph.
- *
- * Each kmer suffix shows up as a distinct SeqVertex, attached in the same structure as in the kmer
- * graph. Nodes that are sources are mapped to SeqVertex nodes that contain all of their sequence
- *
- * @return a newly allocated SequenceGraph
- */
- // TODO -- should override base class method
+ @Override
public SeqGraph convertToSequenceGraph() {
buildGraphIfNecessary();
-
- final SeqGraph seqGraph = new SeqGraph(kmerSize);
- final Map vertexMap = new HashMap<>();
-
-
- // create all of the equivalent seq graph vertices
- for ( final MultiDeBruijnVertex dv : vertexSet() ) {
- final SeqVertex sv = new SeqVertex(dv.getAdditionalSequence(isSource(dv)));
- sv.setAdditionalInfo(dv.additionalInfo());
- vertexMap.put(dv, sv);
- seqGraph.addVertex(sv);
- }
-
- // walk through the nodes and connect them to their equivalent seq vertices
- for( final MultiSampleEdge e : edgeSet() ) {
- final SeqVertex seqInV = vertexMap.get(getEdgeSource(e));
- final SeqVertex seqOutV = vertexMap.get(getEdgeTarget(e));
- //logger.info("Adding edge " + seqInV + " -> " + seqOutV);
- seqGraph.addEdge(seqInV, seqOutV, new BaseEdge(e.isRef(), e.getMultiplicity()));
- }
-
- return seqGraph;
+ return super.convertToSequenceGraph();
}
private void increaseCountsInMatchedKmers(final SequenceForKmers seqForKmers,
@@ -749,15 +711,15 @@ public class ReadThreadingGraph extends DanglingChainMergingGraph implements Kme
}
private static String pathElementId(final String element) {
- final int parentesysPos = element.indexOf('(');
+ final int openBracketPosition = element.indexOf('(');
- if (parentesysPos == -1)
+ if (openBracketPosition == -1)
return null;
- final int closeParentesysPos = element.lastIndexOf(')');
- if (closeParentesysPos == -1)
+ final int closeBracketPosition = element.lastIndexOf(')');
+ if (closeBracketPosition == -1)
throw new IllegalArgumentException("non-closed id parantesys found in element: " + element);
- final String result = element.substring(parentesysPos + 1,closeParentesysPos).trim();
+ final String result = element.substring(openBracketPosition + 1,closeBracketPosition).trim();
if (result.isEmpty())
throw new IllegalArgumentException("empty id found in element: " + element);
return result;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java
index 1f355359d..f16399e62 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java
@@ -247,7 +247,7 @@ public class VariantDataManager {
logger.warn( "WARNING: Training with very few variant sites! Please check the model reporting PDF to ensure the quality of the model is reliable." );
} else if( trainingData.size() > VRAC.MAX_NUM_TRAINING_DATA ) {
logger.warn( "WARNING: Very large training set detected. Downsampling to " + VRAC.MAX_NUM_TRAINING_DATA + " training variants." );
- Collections.shuffle(trainingData);
+ Collections.shuffle(trainingData, GenomeAnalysisEngine.getRandomGenerator());
return trainingData.subList(0, VRAC.MAX_NUM_TRAINING_DATA);
}
return trainingData;
@@ -295,13 +295,13 @@ public class VariantDataManager {
public List getRandomDataForPlotting( final int numToAdd, final List trainingData, final List antiTrainingData, final List evaluationData ) {
final List returnData = new ExpandingArrayList<>();
- Collections.shuffle(trainingData);
- Collections.shuffle(antiTrainingData);
- Collections.shuffle(evaluationData);
+ Collections.shuffle(trainingData, GenomeAnalysisEngine.getRandomGenerator());
+ Collections.shuffle(antiTrainingData, GenomeAnalysisEngine.getRandomGenerator());
+ Collections.shuffle(evaluationData, GenomeAnalysisEngine.getRandomGenerator());
returnData.addAll(trainingData.subList(0, Math.min(numToAdd, trainingData.size())));
returnData.addAll(antiTrainingData.subList(0, Math.min(numToAdd, antiTrainingData.size())));
returnData.addAll(evaluationData.subList(0, Math.min(numToAdd, evaluationData.size())));
- Collections.shuffle(returnData);
+ Collections.shuffle(returnData, GenomeAnalysisEngine.getRandomGenerator());
return returnData;
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CalculateGenotypePosteriors.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CalculateGenotypePosteriors.java
index 3d1a9da57..b8a0585a2 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CalculateGenotypePosteriors.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CalculateGenotypePosteriors.java
@@ -50,6 +50,7 @@ import org.broadinstitute.sting.commandline.ArgumentCollection;
import org.broadinstitute.sting.commandline.Input;
import org.broadinstitute.sting.commandline.Output;
import org.broadinstitute.sting.commandline.RodBinding;
+import org.broadinstitute.sting.gatk.CommandLineGATK;
import org.broadinstitute.sting.gatk.arguments.StandardVariantContextInputArgumentCollection;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
@@ -59,6 +60,8 @@ import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine;
import org.broadinstitute.sting.gatk.walkers.haplotypecaller.PairHMMLikelihoodCalculationEngine;
import org.broadinstitute.sting.utils.SampleUtils;
import org.broadinstitute.sting.utils.exceptions.UserException;
+import org.broadinstitute.sting.utils.help.DocumentedGATKFeature;
+import org.broadinstitute.sting.utils.help.HelpConstants;
import org.broadinstitute.sting.utils.variant.GATKVCFUtils;
import org.broadinstitute.sting.utils.variant.GATKVariantContextUtils;
import org.broadinstitute.variant.variantcontext.VariantContext;
@@ -148,6 +151,7 @@ import java.util.*;
*
*
*/
+@DocumentedGATKFeature( groupName = HelpConstants.DOCS_CAT_VARMANIP, extraDocs = {CommandLineGATK.class} )
public class CalculateGenotypePosteriors extends RodWalker {
/**
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFs.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFs.java
index 0f577cb23..580380513 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFs.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFs.java
@@ -134,6 +134,9 @@ public class CombineGVCFs extends RodWalker 1 ? startingStates.refBases[1] : (byte)'N');
}
@@ -289,7 +292,8 @@ public class CombineGVCFs extends RodWalker attrs = new HashMap<>(1);
- attrs.put(VCFConstants.END_KEY, Integer.toString(end));
+ if ( !USE_BP_RESOLUTION )
+ attrs.put(VCFConstants.END_KEY, Integer.toString(end));
// genotypes
final GenotypesContext genotypes = GenotypesContext.create();
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java
new file mode 100644
index 000000000..ea93ebe4a
--- /dev/null
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java
@@ -0,0 +1,517 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012 Broad Institute, Inc.
+* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 4. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 5. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 6. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 7. MISCELLANEOUS
+* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.sting.utils.pairhmm;
+
+import com.google.java.contract.Ensures;
+import com.google.java.contract.Requires;
+import org.broadinstitute.sting.utils.QualityUtils;
+
+import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
+import org.broadinstitute.sting.utils.haplotype.Haplotype;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
+import org.broadinstitute.variant.variantcontext.Allele;
+import org.broadinstitute.sting.utils.exceptions.UserException;
+import static org.broadinstitute.sting.utils.pairhmm.PairHMMModel.*;
+
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.BufferedWriter;
+import java.util.Map;
+import java.util.HashMap;
+import java.io.IOException;
+
+
+/**
+ * Created with IntelliJ IDEA.
+ * User: rpoplin, carneiro
+ * Date: 10/16/12
+ */
+public class DebugJNILoglessPairHMM extends LoglessPairHMM {
+
+ private static final boolean dumpSandboxOnly = false; //simulates ifdef
+ private static final boolean debug = false; //simulates ifdef
+ private static final boolean verify = !dumpSandboxOnly && (debug || true); //simulates ifdef
+ private static final boolean debug0_1 = false; //simulates ifdef
+ private static final boolean debug1 = false; //simulates ifdef
+ private static final boolean debug2 = false;
+ private static final boolean debug3 = false;
+
+ //Debugging stats
+ private int numCalls = 0;
+ private int numComputeLikelihoodCalls = 0;
+ protected HashMap filenameToWriter = new HashMap();
+
+ private JNILoglessPairHMM jniPairHMM = null;
+ public DebugJNILoglessPairHMM(final PairHMM.HMM_IMPLEMENTATION hmmType) {
+ super();
+ switch(hmmType) {
+ case VECTOR_LOGLESS_CACHING:
+ jniPairHMM = new VectorLoglessPairHMM();
+ break;
+ default:
+ throw new UserException.BadArgumentValue("pairHMM","Specified JNIPairHMM implementation is unrecognized or incompatible with the HaplotypeCaller. Acceptable options are VECTOR_LOGLESS_CACHING");
+ }
+ }
+
+ @Override
+ public void close()
+ {
+ jniPairHMM.close();
+ debugClose();
+ }
+
+ //Used only when testing parts of the compute kernel
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public void initialize( final int readMaxLength, final int haplotypeMaxLength ) {
+ if(verify)
+ super.initialize(readMaxLength, haplotypeMaxLength);
+ if(debug3)
+ {
+ System.out.println("Java: alloc initialized readMaxLength : "+readMaxLength+" haplotypeMaxLength : "+haplotypeMaxLength);
+ debugDump("lengths_java.txt", String.format("%d %d\n",readMaxLength, haplotypeMaxLength),
+ true);
+ }
+ if(debug2)
+ jniInitialize(readMaxLength, haplotypeMaxLength);
+ }
+
+ private HashMap haplotypeToHaplotypeListIdxMap = null;
+ //Used to transfer data to JNI
+ //Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public void initialize( final List haplotypes, final Map> perSampleReadList,
+ final int readMaxLength, final int haplotypeMaxLength ) {
+ if(verify)
+ {
+ super.initialize(haplotypes, perSampleReadList, readMaxLength, haplotypeMaxLength);
+ jniPairHMM.initialize(haplotypes, perSampleReadList, readMaxLength, haplotypeMaxLength);
+ haplotypeToHaplotypeListIdxMap = jniPairHMM.getHaplotypeToHaplotypeListIdxMap();
+ }
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public void finalizeRegion()
+ {
+ if(!dumpSandboxOnly)
+ jniPairHMM.finalizeRegion();
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public PerReadAlleleLikelihoodMap computeLikelihoods( final List reads, final Map alleleHaplotypeMap, final Map GCPArrayMap ) {
+ // (re)initialize the pairHMM only if necessary
+ final int readMaxLength = verify ? findMaxReadLength(reads) : 0;
+ final int haplotypeMaxLength = verify ? findMaxHaplotypeLength(alleleHaplotypeMap) : 0;
+ if(verify)
+ {
+ if (!initialized || readMaxLength > maxReadLength || haplotypeMaxLength > maxHaplotypeLength)
+ { initialize(readMaxLength, haplotypeMaxLength); }
+ if ( ! initialized )
+ throw new IllegalStateException("Must call initialize before calling jniComputeLikelihoods in debug/verify mode");
+ }
+ int readListSize = reads.size();
+ int numHaplotypes = alleleHaplotypeMap.size();
+ int numTestcases = readListSize*numHaplotypes;
+ if(debug0_1)
+ System.out.println("Java numReads "+readListSize+" numHaplotypes "+numHaplotypes);
+ int idx = 0;
+ for(GATKSAMRecord read : reads)
+ {
+ byte [] overallGCP = GCPArrayMap.get(read);
+ if(debug0_1)
+ System.out.println("Java read length "+read.getReadBases().length);
+ if(debug3)
+ {
+ for(int i=0;i currEntry : alleleHaplotypeMap.entrySet()) //order is important - access in same order always
+ {
+ byte[] haplotypeBases = currEntry.getValue().getBases();
+ if(debug0_1)
+ System.out.println("Java haplotype length "+haplotypeBases.length);
+ if(debug3)
+ {
+ for(int i=0;i currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always
+ {
+ idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(currEntry.getValue());
+ likelihoodArray[idx] = tmpArray[idxInsideHaplotypeList];
+ ++idx;
+ }
+ readIdx += numHaplotypes;
+ }
+ //for floating point values, no exact equality
+ //check whether numbers are close in terms of abs_error or relative_error
+ //For very large values, relative_error is relevant
+ //For very small values, abs_error is relevant
+ for(int i=0;i 1e-5 && relative_error > 1e-5)
+ {
+ toDump = true;
+ break;
+ }
+ }
+ }
+ //if numbers are not close, then dump out the data that produced the inconsistency
+ if(toDump)
+ {
+ idx = 0;
+ System.out.println("Dump : Java numReads "+readListSize+" numHaplotypes "+numHaplotypes);
+ boolean firstLine = true;
+ for(GATKSAMRecord read : reads)
+ {
+ byte [] overallGCP = GCPArrayMap.get(read);
+ byte[] tmpByteArray = new byte[read.getReadBases().length];
+ for (Map.Entry currEntry : alleleHaplotypeMap.entrySet()) //order is important - access in same order always
+ {
+ byte[] haplotypeBases = currEntry.getValue().getBases();
+ debugDump("debug_dump.txt",new String(haplotypeBases)+" ",true);
+ debugDump("debug_dump.txt",new String(read.getReadBases())+" ",true);
+ for(int k=0;k currEntry : filenameToWriter.entrySet()) {
+ BufferedWriter currWriter = currEntry.getValue();
+ try
+ {
+ currWriter.flush();
+ currWriter.close();
+ }
+ catch(IOException e)
+ {
+ e.printStackTrace();
+
+ }
+ }
+ filenameToWriter.clear();
+ }
+}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java
new file mode 100644
index 000000000..f039cc295
--- /dev/null
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java
@@ -0,0 +1,63 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012 Broad Institute, Inc.
+* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 4. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 5. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 6. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 7. MISCELLANEOUS
+* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.sting.utils.pairhmm;
+
+import org.broadinstitute.sting.utils.haplotype.Haplotype;
+
+import java.util.HashMap;
+
+
+/**
+ * Created with IntelliJ IDEA.
+ * User: rpoplin, carneiro
+ * Date: 10/16/12
+ */
+public abstract class JNILoglessPairHMM extends LoglessPairHMM {
+ public abstract HashMap getHaplotypeToHaplotypeListIdxMap();
+ protected long setupTime = 0;
+
+}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/LoglessPairHMM.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/LoglessPairHMM.java
index ed35e6970..31a0d1363 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/LoglessPairHMM.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/LoglessPairHMM.java
@@ -66,7 +66,7 @@ public class LoglessPairHMM extends N2MemoryPairHMM {
protected static final double TRISTATE_CORRECTION = 3.0;
-
+
/**
* {@inheritDoc}
*/
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java
new file mode 100644
index 000000000..e69d9ea50
--- /dev/null
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java
@@ -0,0 +1,340 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012 Broad Institute, Inc.
+* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 4. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 5. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 6. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 7. MISCELLANEOUS
+* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.sting.utils.pairhmm;
+
+import com.google.java.contract.Ensures;
+import com.google.java.contract.Requires;
+import org.broadinstitute.sting.utils.QualityUtils;
+
+import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
+import org.broadinstitute.sting.utils.haplotype.Haplotype;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
+import org.broadinstitute.variant.variantcontext.Allele;
+
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+
+//For loading library from jar
+import java.io.File;
+import java.io.FileNotFoundException;
+import java.io.FileOutputStream;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.OutputStream;
+
+/**
+ * Created with IntelliJ IDEA.
+ * User: rpoplin, carneiro
+ * Date: 10/16/12
+ */
+public class VectorLoglessPairHMM extends JNILoglessPairHMM {
+
+ //For machine capabilities
+ public static final long sse41Mask = 1;
+ public static final long sse42Mask = 2;
+ public static final long avxMask = 4;
+ public static final long enableAll = 0xFFFFFFFFFFFFFFFFl;
+
+ //Used to copy references to byteArrays to JNI from reads
+ protected class JNIReadDataHolderClass {
+ public byte[] readBases = null;
+ public byte[] readQuals = null;
+ public byte[] insertionGOP = null;
+ public byte[] deletionGOP = null;
+ public byte[] overallGCP = null;
+ }
+
+ //Used to copy references to byteArrays to JNI from haplotypes
+ protected class JNIHaplotypeDataHolderClass {
+ public byte[] haplotypeBases = null;
+ }
+
+ /**
+ * Return 64-bit mask representing machine capabilities
+ * Bit 0 is LSB, bit 63 MSB
+ * Bit 0 represents sse4.1 availability
+ * Bit 1 represents sse4.2 availability
+ * Bit 2 represents AVX availability
+ */
+ public native long jniGetMachineType();
+
+ /**
+ * Function to initialize the fields of JNIReadDataHolderClass and JNIHaplotypeDataHolderClass from JVM.
+ * C++ codegets FieldIDs for these classes once and re-uses these IDs for the remainder of the program. Field IDs do not
+ * change per JVM session
+ * @param readDataHolderClass class type of JNIReadDataHolderClass
+ * @param haplotypeDataHolderClass class type of JNIHaplotypeDataHolderClass
+ * @param mask mask is a 64 bit integer identical to the one received from jniGetMachineType(). Users can disable usage of some hardware features by zeroing some bits in the mask
+ * */
+ private native void jniInitializeClassFieldsAndMachineMask(Class readDataHolderClass, Class haplotypeDataHolderClass, long mask);
+
+ private static Boolean isVectorLoglessPairHMMLibraryLoaded = false;
+ //The constructor is called only once inside PairHMMLikelihoodCalculationEngine
+ public VectorLoglessPairHMM() {
+ super();
+
+ logger.warn("WARNING: the VectorLoglessPairHMM is an experimental implementation still under active development. " +
+ "Use at your own risk!");
+
+ synchronized(isVectorLoglessPairHMMLibraryLoaded) {
+ //Load the library and initialize the FieldIDs
+ if(!isVectorLoglessPairHMMLibraryLoaded) {
+ try
+ {
+ //Try loading from Java's library path first
+ //Useful if someone builds his/her own library and wants to override the bundled
+ //implementation without modifying the Java code
+ System.loadLibrary("VectorLoglessPairHMM");
+ logger.info("libVectorLoglessPairHMM found in JVM library path");
+ }
+ catch(UnsatisfiedLinkError ule)
+ {
+ //Could not load from Java's library path - try unpacking from jar
+ try
+ {
+ logger.debug("libVectorLoglessPairHMM not found in JVM library path - trying to unpack from StingUtils.jar");
+ loadLibraryFromJar("/org/broadinstitute/sting/utils/pairhmm/libVectorLoglessPairHMM.so");
+ logger.debug("libVectorLoglessPairHMM unpacked successfully from StingUtils.jar");
+ }
+ catch(IOException ioe)
+ {
+ //Throw the UnsatisfiedLinkError to make it clear to the user what failed
+ throw ule;
+ }
+ }
+
+ isVectorLoglessPairHMMLibraryLoaded = true;
+ jniInitializeClassFieldsAndMachineMask(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class, enableAll); //need to do this only once
+ }
+ }
+ }
+
+ private native void jniInitializeHaplotypes(final int numHaplotypes, JNIHaplotypeDataHolderClass[] haplotypeDataArray);
+ //Hold the mapping between haplotype and index in the list of Haplotypes passed to initialize
+ //Use this mapping in computeLikelihoods to find the likelihood value corresponding to a given Haplotype
+ private HashMap haplotypeToHaplotypeListIdxMap = new HashMap();
+ private JNIHaplotypeDataHolderClass[] mHaplotypeDataArray;
+ @Override
+ public HashMap getHaplotypeToHaplotypeListIdxMap() { return haplotypeToHaplotypeListIdxMap; }
+
+ //Used to transfer data to JNI
+ //Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public void initialize( final List haplotypes, final Map> perSampleReadList,
+ final int readMaxLength, final int haplotypeMaxLength ) {
+ int numHaplotypes = haplotypes.size();
+ mHaplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes];
+ int idx = 0;
+ haplotypeToHaplotypeListIdxMap.clear();
+ for(final Haplotype currHaplotype : haplotypes)
+ {
+ mHaplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass();
+ mHaplotypeDataArray[idx].haplotypeBases = currHaplotype.getBases();
+ haplotypeToHaplotypeListIdxMap.put(currHaplotype, idx);
+ ++idx;
+ }
+ jniInitializeHaplotypes(numHaplotypes, mHaplotypeDataArray);
+ }
+ /**
+ * Tell JNI to release arrays - really important if native code is directly accessing Java memory, if not
+ * accessing Java memory directly, still important to release memory from C++
+ */
+ private native void jniFinalizeRegion();
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public void finalizeRegion()
+ {
+ jniFinalizeRegion();
+ }
+
+ /**
+ * Real compute kernel
+ */
+ private native void jniComputeLikelihoods(int numReads, int numHaplotypes, JNIReadDataHolderClass[] readDataArray,
+ JNIHaplotypeDataHolderClass[] haplotypeDataArray, double[] likelihoodArray, int maxNumThreadsToUse);
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public PerReadAlleleLikelihoodMap computeLikelihoods( final List reads, final Map alleleHaplotypeMap, final Map GCPArrayMap ) {
+ if(doProfiling)
+ startTime = System.nanoTime();
+ int readListSize = reads.size();
+ int numHaplotypes = alleleHaplotypeMap.size();
+ int numTestcases = readListSize*numHaplotypes;
+ JNIReadDataHolderClass[] readDataArray = new JNIReadDataHolderClass[readListSize];
+ int idx = 0;
+ for(GATKSAMRecord read : reads)
+ {
+ readDataArray[idx] = new JNIReadDataHolderClass();
+ readDataArray[idx].readBases = read.getReadBases();
+ readDataArray[idx].readQuals = read.getBaseQualities();
+ readDataArray[idx].insertionGOP = read.getBaseInsertionQualities();
+ readDataArray[idx].deletionGOP = read.getBaseDeletionQualities();
+ readDataArray[idx].overallGCP = GCPArrayMap.get(read);
+ ++idx;
+ }
+
+ mLikelihoodArray = new double[readListSize*numHaplotypes]; //to store results
+ if(doProfiling)
+ setupTime += (System.nanoTime() - startTime);
+ //for(reads)
+ // for(haplotypes)
+ // compute_full_prob()
+ jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, mHaplotypeDataArray, mLikelihoodArray, 12);
+
+ final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap();
+ idx = 0;
+ int idxInsideHaplotypeList = 0;
+ int readIdx = 0;
+ for(GATKSAMRecord read : reads)
+ {
+ for (Map.Entry currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always
+ {
+ //Since the order of haplotypes in the List and alleleHaplotypeMap is different,
+ //get idx of current haplotype in the list and use this idx to get the right likelihoodValue
+ idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(currEntry.getValue());
+ likelihoodMap.add(read, currEntry.getKey(), mLikelihoodArray[readIdx + idxInsideHaplotypeList]);
+ ++idx;
+ }
+ readIdx += numHaplotypes;
+ }
+ if(doProfiling)
+ computeTime += (System.nanoTime() - startTime);
+ return likelihoodMap;
+ }
+
+ /**
+ * Print final profiling information from native code
+ */
+ public native void jniClose();
+ @Override
+ public void close()
+ {
+ System.out.println("Time spent in setup for JNI call : "+(setupTime*1e-9));
+ super.close();
+ jniClose();
+ }
+
+ //Copied from http://frommyplayground.com/how-to-load-native-jni-library-from-jar
+ /**
+ * Loads library from current JAR archive
+ *
+ * The file from JAR is copied into system temporary directory and then loaded. The temporary file is deleted after exiting.
+ * Method uses String as filename because the pathname is "abstract", not system-dependent.
+ *
+ * @param path The filename inside JAR as absolute path (beginning with '/'), e.g. /package/File.ext
+ * @throws IOException If temporary file creation or read/write operation fails
+ * @throws IllegalArgumentException If source file (param path) does not exist
+ * @throws IllegalArgumentException If the path is not absolute or if the filename is shorter than three characters (restriction of {@see File#createTempFile(java.lang.String, java.lang.String)}).
+ */
+ public static void loadLibraryFromJar(String path) throws IOException {
+
+ if (!path.startsWith("/")) {
+ throw new IllegalArgumentException("The path to be absolute (start with '/').");
+ }
+
+ // Obtain filename from path
+ String[] parts = path.split("/");
+ String filename = (parts.length > 1) ? parts[parts.length - 1] : null;
+
+ // Split filename to prexif and suffix (extension)
+ String prefix = "";
+ String suffix = null;
+ if (filename != null) {
+ parts = filename.split("\\.", 2);
+ prefix = parts[0];
+ suffix = (parts.length > 1) ? "."+parts[parts.length - 1] : null; // Thanks, davs! :-)
+ }
+
+ // Check if the filename is okay
+ if (filename == null || prefix.length() < 3) {
+ throw new IllegalArgumentException("The filename has to be at least 3 characters long.");
+ }
+
+ // Prepare temporary file
+ File temp = File.createTempFile(prefix, suffix);
+ //System.out.println("Temp lib file "+temp.getAbsolutePath());
+ temp.deleteOnExit();
+
+ if (!temp.exists()) {
+ throw new FileNotFoundException("File " + temp.getAbsolutePath() + " does not exist.");
+ }
+
+ // Prepare buffer for data copying
+ byte[] buffer = new byte[1024];
+ int readBytes;
+
+ // Open and check input stream
+ InputStream is = VectorLoglessPairHMM.class.getResourceAsStream(path);
+ if (is == null) {
+ throw new FileNotFoundException("File " + path + " was not found inside JAR.");
+ }
+
+ // Open output stream and copy data between source file in JAR and the temporary file
+ OutputStream os = new FileOutputStream(temp);
+ try {
+ while ((readBytes = is.read(buffer)) != -1) {
+ os.write(buffer, 0, readBytes);
+ }
+ } finally {
+ // If read/write fails, close streams safely before throwing an exception
+ os.close();
+ is.close();
+ }
+
+ // Finally, load the library
+ System.load(temp.getAbsolutePath());
+ }
+}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java
index b1c280748..0ec1dd996 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java
@@ -46,6 +46,7 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
+import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.utils.MannWhitneyU;
import org.testng.Assert;
import org.testng.annotations.BeforeClass;
@@ -75,9 +76,9 @@ public class RankSumUnitTest {
makeDistribution(distribution20_40, 40, skew, observations/2);
// shuffle the observations
- Collections.shuffle(distribution20);
- Collections.shuffle(distribution30);
- Collections.shuffle(distribution20_40);
+ Collections.shuffle(distribution20, GenomeAnalysisEngine.getRandomGenerator());
+ Collections.shuffle(distribution30, GenomeAnalysisEngine.getRandomGenerator());
+ Collections.shuffle(distribution20_40, GenomeAnalysisEngine.getRandomGenerator());
}
private static void makeDistribution(final List result, final int target, final int skew, final int numObservations) {
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/StrandOddsRatioUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/StrandOddsRatioUnitTest.java
new file mode 100644
index 000000000..562736f0a
--- /dev/null
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/StrandOddsRatioUnitTest.java
@@ -0,0 +1,100 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012 Broad Institute, Inc.
+* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 4. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 5. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 6. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 7. MISCELLANEOUS
+* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.sting.gatk.walkers.annotator;
+
+import org.testng.Assert;
+import org.testng.annotations.DataProvider;
+import org.testng.annotations.Test;
+
+import java.util.ArrayList;
+import java.lang.Integer;
+import java.util.List;
+
+/**
+ * Created by haasb on 3/5/14.
+ */
+public class StrandOddsRatioUnitTest {
+ private static double DELTA_PRECISION = 0.001;
+
+ @DataProvider(name = "UsingSOR")
+ public Object[][] makeUsingSORData() {
+ List tests = new ArrayList<>();
+ tests.add(new Object[]{0, 0, 0, 0, 2.0});
+ tests.add(new Object[]{100000, 100000, 100000, 100000, 2.0} );
+ tests.add(new Object[]{Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 2.0} );
+
+ tests.add(new Object[]{0, 0, 100000, 100000, 2.0});
+ tests.add(new Object[]{0, 0, Integer.MAX_VALUE, Integer.MAX_VALUE, 2.0});
+
+ tests.add(new Object[]{100000,100000,100000,0, 1000001.000001});
+ tests.add(new Object[]{100,100,100,0, 1001.000999});
+ tests.add(new Object[]{Integer.MAX_VALUE,Integer.MAX_VALUE,Integer.MAX_VALUE,0, 21474836471.0});
+
+ tests.add(new Object[]{13736,9047,41,1433, 52.95947});
+ tests.add(new Object[]{66, 14, 64, 4, 3.63482});
+ tests.add(new Object[]{351169, 306836, 153739, 2379, 56.48043});
+ tests.add(new Object[]{116449, 131216, 289, 16957, 52.07302});
+ tests.add(new Object[]{137, 159, 9, 23, 2.64460});
+ tests.add(new Object[]{129, 90, 21, 20, 2.09757});
+ tests.add(new Object[]{14054, 9160, 16, 7827, 745.89657});
+ tests.add(new Object[]{32803, 9184, 32117, 3283, 3.10399});
+ tests.add(new Object[]{2068, 6796, 1133, 0, 37235.43791});
+
+ return tests.toArray(new Object[][]{});
+ }
+
+ @Test(dataProvider = "UsingSOR")
+ public void testUsingSOR(final int refpos, final int refneg, final int altpos, final int altneg, double expectedOddsRatio ) {
+ int[][] contingencyTable = new int[2][2];
+ contingencyTable[0][0] = refpos;
+ contingencyTable[0][1] = refneg;
+ contingencyTable[1][0] = altpos;
+ contingencyTable[1][1] = altneg;
+ final double ratio = new StrandOddsRatio().symmetricOddsRatio(contingencyTable);
+ Assert.assertEquals(ratio, expectedOddsRatio, DELTA_PRECISION, "Pass");
+ }
+}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceIntegrationTest.java
index a650e0f6f..58ff0a7b3 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceIntegrationTest.java
@@ -47,38 +47,63 @@
package org.broadinstitute.sting.gatk.walkers.fasta;
import org.broadinstitute.sting.WalkerTest;
+import org.broadinstitute.sting.utils.exceptions.UserException;
import org.testng.annotations.Test;
import java.util.Arrays;
public class FastaAlternateReferenceIntegrationTest extends WalkerTest {
+
@Test
- public void testIntervals() {
+ public void testReferenceOnly() {
- String md5_1 = "328d2d52cedfdc52da7d1abff487633d";
-
- WalkerTestSpec spec1a = new WalkerTestSpec(
- "-T FastaAlternateReferenceMaker -R " + b36KGReference + " -L 1:10,000,100-10,000,500 -L 1:10,100,000-10,101,000 -L 1:10,900,000-10,900,001 -o %s",
- 1,
- Arrays.asList(md5_1));
- executeTest("testFastaReference", spec1a);
-
- WalkerTestSpec spec1b = new WalkerTestSpec(
+ WalkerTestSpec spec = new WalkerTestSpec(
"-T FastaReferenceMaker -R " + b36KGReference + " -L 1:10,000,100-10,000,500 -L 1:10,100,000-10,101,000 -L 1:10,900,000-10,900,001 -o %s",
1,
- Arrays.asList(md5_1));
- executeTest("testFastaReference", spec1b);
+ Arrays.asList("328d2d52cedfdc52da7d1abff487633d"));
+ executeTest("test FastaReference", spec);
+ }
- WalkerTestSpec spec2 = new WalkerTestSpec(
+ @Test
+ public void testIndelsAndSnpMask() {
+
+ WalkerTestSpec spec = new WalkerTestSpec(
"-T FastaAlternateReferenceMaker -R " + b36KGReference + " -V " + validationDataLocation + "NA12878.chr1_10mb_11mb.slx.indels.vcf4 --snpmask:vcf " + b36dbSNP129 + " -L 1:10,075,000-10,075,380 -L 1:10,093,447-10,093,847 -L 1:10,271,252-10,271,452 -o %s",
1,
Arrays.asList("ef481be9962e21d09847b8a1d4a4ff65"));
- executeTest("testFastaAlternateReferenceIndels", spec2);
+ executeTest("test indels", spec);
+ }
- WalkerTestSpec spec3 = new WalkerTestSpec(
+ @Test
+ public void testSnps() {
+
+ WalkerTestSpec spec = new WalkerTestSpec(
"-T FastaAlternateReferenceMaker -R " + b36KGReference + " -V " + GATKDataLocation + "dbsnp_129_b36.vcf -L 1:10,023,400-10,023,500 -L 1:10,029,200-10,029,500 -o %s",
1,
Arrays.asList("8b6cd2e20c381f9819aab2d270f5e641"));
- executeTest("testFastaAlternateReferenceSnps", spec3);
+ executeTest("test SNPs", spec);
+ }
+
+ @Test
+ public void testBadIupacInput() {
+
+ // cannot use 'expectedExceptions = UserException.BadInput.class' because it technically gets thrown as a RuntimeException by the engine
+ try {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ "-T FastaAlternateReferenceMaker -R " + b36KGReference + " --useIUPAC -V " + GATKDataLocation + "dbsnp_129_b36.vcf -L 1:10,023,400-10,023,500 -L 1:10,029,200-10,029,500 -o %s",
+ 1,
+ Arrays.asList("FAILFAILFAILFAILFAILFAILFAILFAIL"));
+ executeTest("test bad input", spec);
+ } catch (Exception e) {} // do nothing
+ }
+
+ @Test
+ public void testIupac() {
+
+ WalkerTestSpec spec = new WalkerTestSpec(
+ "-T FastaAlternateReferenceMaker -R " + b37KGReference + " --useIUPAC -V " + privateTestDir + "NA12878.WGS.b37.chr20.firstMB.vcf -L 20:61050-66380 -o %s",
+ 1,
+ Arrays.asList("5feb2a576ff2ed1745a007eaa36448b3"));
+ executeTest("test iupac", spec);
}
}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerComplexAndSymbolicVariantsIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerComplexAndSymbolicVariantsIntegrationTest.java
index 2838648d5..24e47879d 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerComplexAndSymbolicVariantsIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerComplexAndSymbolicVariantsIntegrationTest.java
@@ -94,6 +94,19 @@ public class HaplotypeCallerComplexAndSymbolicVariantsIntegrationTest extends Wa
@Test
public void testHaplotypeCallerMultiSampleGGAMultiAllelic() {
HCTestComplexGGA(NA12878_CHR20_BAM, "-L 20:133041-133161 -L 20:300207-300337",
- "f50e0b35e2240b19b1b8b6dfa0cf9796");
+ "5ac3bfe1da1d411b52a98ef3debbd318");
}
+
+ private void HCTestComplexConsensusMode(String bam, String args, String md5) {
+ final String base = String.format("-T HaplotypeCaller --disableDithering --pcr_indel_model NONE -R %s -I %s", REF, bam) + " --no_cmdline_in_header -o %s -consensus -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf -alleles " + validationDataLocation + "phase1.projectConsensus.chr20.raw.snps.vcf";
+ final WalkerTestSpec spec = new WalkerTestSpec(base + " " + args, Arrays.asList(md5));
+ executeTest("testHaplotypeCallerComplexConsensusMode: args=" + args, spec);
+ }
+
+ @Test
+ public void testHaplotypeCallerMultiSampleConsensusModeComplex() {
+ HCTestComplexGGA(NA12878_CHR20_BAM, "-L 20:119673-119823 -L 20:121408-121538 -L 20:133041-133161 -L 20:300207-300337",
+ "61972c7c0d378e756f3b4d99aed9d0cf");
+ }
+
}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java
index 8ca67f31d..bde1e7ef5 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java
@@ -68,8 +68,8 @@ public class HaplotypeCallerGVCFIntegrationTest extends WalkerTest {
// this functionality can be adapted to provide input data for whatever you might want in your data
tests.add(new Object[]{NA12878_PCRFREE, HaplotypeCaller.ReferenceConfidenceMode.NONE, PCRFreeIntervals, "50323a284788c8220c9226037c7003b5"});
- tests.add(new Object[]{NA12878_PCRFREE, HaplotypeCaller.ReferenceConfidenceMode.BP_RESOLUTION, PCRFreeIntervals, "7c16aa8e35de9f418533efac3bae6551"});
- tests.add(new Object[]{NA12878_PCRFREE, HaplotypeCaller.ReferenceConfidenceMode.GVCF, PCRFreeIntervals, "7e1e193d70187774f9740d475e0f1cc1"});
+ tests.add(new Object[]{NA12878_PCRFREE, HaplotypeCaller.ReferenceConfidenceMode.BP_RESOLUTION, PCRFreeIntervals, "96fea2caf0a40df3feb268e8b14da670"});
+ tests.add(new Object[]{NA12878_PCRFREE, HaplotypeCaller.ReferenceConfidenceMode.GVCF, PCRFreeIntervals, "19efc8020f31d1b68d80c50df0629e50"});
tests.add(new Object[]{NA12878_WEx, HaplotypeCaller.ReferenceConfidenceMode.NONE, WExIntervals, "39bf5fe3911d0c646eefa8f79894f4df"});
tests.add(new Object[]{NA12878_WEx, HaplotypeCaller.ReferenceConfidenceMode.BP_RESOLUTION, WExIntervals, "d926d653500a970280ad7828d9ee2b84"});
tests.add(new Object[]{NA12878_WEx, HaplotypeCaller.ReferenceConfidenceMode.GVCF, WExIntervals, "83ddc16e4f0900429b2da30e582994aa"});
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
index 386fc3800..8db4d5066 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
@@ -227,7 +227,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
public void HCTestDBSNPAnnotationWGS() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T HaplotypeCaller --disableDithering --pcr_indel_model NONE -R " + b37KGReference + " --no_cmdline_in_header -I " + NA12878_PCRFREE + " -o %s -L 20:10,000,000-10,100,000 -D " + b37dbSNP132, 1,
- Arrays.asList("0998be22d7af4372247f5a0338f9446b"));
+ Arrays.asList("7c3254ead383e2b9a51b242f6de2a5b2"));
executeTest("HC calling with dbSNP ID annotation on WGS intervals", spec);
}
@@ -244,7 +244,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
public void HCTestDBSNPAnnotationWGSGraphBased() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T HaplotypeCaller -likelihoodEngine GraphBased --disableDithering --pcr_indel_model NONE -R " + b37KGReference + " --no_cmdline_in_header -I " + NA12878_PCRFREE + " -o %s -L 20:10,000,000-10,100,000 -D " + b37dbSNP132, 1,
- Arrays.asList("1aeed297a3cb41940d83eac499a2ce07"));
+ Arrays.asList("eda8f91091fe462205d687ec49fc61e7"));
executeTest("HC calling with dbSNP ID annotation on WGS intervals", spec);
}
@@ -276,7 +276,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
public void HCTestAggressivePcrIndelModelWGS() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T HaplotypeCaller --disableDithering --pcr_indel_model AGGRESSIVE -R " + b37KGReference + " --no_cmdline_in_header -I " + NA12878_BAM + " -o %s -L 20:10,000,000-10,300,000", 1,
- Arrays.asList("f426f4c2986e1dea8f3f55951ef8e013"));
+ Arrays.asList("73c52372a1a80f052ea2b728ee17bf22"));
executeTest("HC calling with aggressive indel error modeling on WGS intervals", spec);
}
@@ -284,7 +284,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
public void HCTestConservativePcrIndelModelWGS() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T HaplotypeCaller --disableDithering --pcr_indel_model CONSERVATIVE -R " + b37KGReference + " --no_cmdline_in_header -I " + NA12878_BAM + " -o %s -L 20:10,000,000-10,300,000", 1,
- Arrays.asList("dcb38cb9280f2c3059a09d323db1c633"));
+ Arrays.asList("4e10d49b8af23d5ef3a28cb702d10a4b"));
executeTest("HC calling with conservative indel error modeling on WGS intervals", spec);
}
@@ -298,4 +298,25 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
spec.disableShadowBCF();
executeTest("testGraphBasedNoSuchEdgeBugFix", spec);
}
+
+ @Test
+ public void testLackSensitivityDueToBadHaplotypeSelectionFix() {
+ final String commandLine = String.format("-T HaplotypeCaller -R %s -I %s -L %s --no_cmdline_in_header ",
+ b37KGReferenceWithDecoy, privateTestDir + "hc-lack-sensitivity.bam", privateTestDir + "hc-lack-sensitivity.interval_list",
+ HaplotypeCaller.OPTIMAL_GVCF_INDEX_TYPE, HaplotypeCaller.OPTIMAL_GVCF_INDEX_PARAMETER);
+ final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", Arrays.asList("e2e6647f7c96e91aeead7301017dc800"));
+ spec.disableShadowBCF();
+ executeTest("testLackSensitivityDueToBadHaplotypeSelectionFix", spec);
+ }
+
+ @Test
+ public void testBadLikelihoodsDueToBadHaplotypeSelectionFix() {
+ final String commandLine = String.format("-T HaplotypeCaller -R %s -I %s -L %s --no_cmdline_in_header ",
+ hg19RefereneWithChrPrefixInChromosomeNames, privateTestDir + "bad-likelihoods.bam", privateTestDir + "bad-likelihoods.interval_list",
+ HaplotypeCaller.OPTIMAL_GVCF_INDEX_TYPE, HaplotypeCaller.OPTIMAL_GVCF_INDEX_PARAMETER);
+ final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", Arrays.asList("cbda30145523bf05e0413157f1a00b3e"));
+ spec.disableShadowBCF();
+ executeTest("testBadLikelihoodsDueToBadHaplotypeSelectionFix", spec);
+ }
+
}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerParallelIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerParallelIntegrationTest.java
index 23513f314..b801f05a9 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerParallelIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerParallelIntegrationTest.java
@@ -47,7 +47,6 @@
package org.broadinstitute.sting.gatk.walkers.haplotypecaller;
import org.broadinstitute.sting.WalkerTest;
-import org.broadinstitute.sting.utils.haplotypeBAMWriter.HaplotypeBAMWriter;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
@@ -61,7 +60,7 @@ public class HaplotypeCallerParallelIntegrationTest extends WalkerTest {
List tests = new ArrayList<>();
for ( final int nct : Arrays.asList(1, 2, 4) ) {
- tests.add(new Object[]{nct, "1f463bf3a06c401006858bc446ecea54"});
+ tests.add(new Object[]{nct, "fd9324a574f9204f7308fc1af422fdcc"});
}
return tests.toArray(new Object[][]{});
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/CommonSuffixMergerUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/CommonSuffixMergerUnitTest.java
index 0ddf7544d..da1474db1 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/CommonSuffixMergerUnitTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/CommonSuffixMergerUnitTest.java
@@ -52,10 +52,11 @@ import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
import java.io.File;
+import java.io.IOException;
import java.util.*;
public class CommonSuffixMergerUnitTest extends BaseTest {
- private final static boolean PRINT_GRAPHS = true;
+ private final static boolean PRINT_GRAPHS = false;
@DataProvider(name = "CompleteCycleData")
public Object[][] makeCompleteCycleData() {
@@ -134,11 +135,35 @@ public class CommonSuffixMergerUnitTest extends BaseTest {
return toUse.toArray(new Object[][]{});
}
+ /**
+ * Compares KBestHaplotype solutions, first by the haplotype base sequence and the by their score.
+ */
+ private static final Comparator KBESTHAPLOTYPE_COMPARATOR = new Comparator() {
+
+ /**
+ * Compares KBestHaplotype solutions, first by the haplotype base sequence and the by their score.
+ *
+ * @return {@inheritDoc}
+ */
+ @Override
+ public int compare(final KBestHaplotype o1,final KBestHaplotype o2) {
+ final int baseCmp = o1.haplotype().getBaseString().compareTo(o2.haplotype().getBaseString());
+ if (baseCmp != 0)
+ return baseCmp;
+ return - Double.compare(o1.score(),o2.score());
+ }
+ };
+
+
public static void assertSameHaplotypes(final String name, final SeqGraph actual, final SeqGraph original) {
+ final KBestHaplotypeFinder originalKBestHaplotypes = new KBestHaplotypeFinder(original,original.getSources(),original.getSinks());
+ final KBestHaplotypeFinder actualKBestHaplotypes = new KBestHaplotypeFinder(actual,actual.getSources(),actual.getSinks());
+ final List sortedOriginalKBestHaplotypes = new ArrayList<>(originalKBestHaplotypes);
+ Collections.sort(sortedOriginalKBestHaplotypes, KBESTHAPLOTYPE_COMPARATOR);
+ final List sortedActualKBestHaplotypes = new ArrayList<>(actualKBestHaplotypes);
+ Collections.sort(sortedActualKBestHaplotypes, KBESTHAPLOTYPE_COMPARATOR);
try {
final Set haplotypes = new HashSet();
- final List originalKBestHaplotypes = new KBestHaplotypeFinder(original,original.getSources(),original.getSinks());
- final List actualKBestHaplotypes = new KBestHaplotypeFinder(actual,actual.getSources(),actual.getSinks());
for (final KBestHaplotype kbh : originalKBestHaplotypes)
haplotypes.add(new String(kbh.bases()));
@@ -148,14 +173,16 @@ public class CommonSuffixMergerUnitTest extends BaseTest {
Assert.assertTrue(haplotypes.contains(h), "Failed to find haplotype " + h);
}
- if ( actualKBestHaplotypes.size() == originalKBestHaplotypes.size() ) {
- for ( int i = 0; i < originalKBestHaplotypes.size(); i++ ) {
- Assert.assertTrue(actualKBestHaplotypes.get(i).haplotype().getBaseString().equals(originalKBestHaplotypes.get(i).haplotype().getBaseString()), "Paths not equal " + actualKBestHaplotypes.get(i).haplotype() + " vs. original " + originalKBestHaplotypes.get(i).haplotype());
- }
- }
+ Assert.assertEquals(sortedActualKBestHaplotypes,sortedOriginalKBestHaplotypes);
} catch ( AssertionError e ) {
if ( PRINT_GRAPHS ) original.printGraph(new File(String.format("%s.original.dot", name, actual.vertexSet().size())), 0);
if ( PRINT_GRAPHS ) actual.printGraph(new File(String.format("%s.actual.dot", name, actual.vertexSet().size())), 0);
+ try {
+ if ( PRINT_GRAPHS ) originalKBestHaplotypes.printDOTFile(String.format("%s.original.finder.dot",name));
+ if ( PRINT_GRAPHS ) actualKBestHaplotypes.printDOTFile(String.format("%s.actual.finder.dot",name));
+ } catch (IOException e2) {
+ // do nothing.
+ }
throw e;
}
}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/KBestHaplotypeFinderUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/KBestHaplotypeFinderUnitTest.java
index 6dc3d5d67..26c511b6e 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/KBestHaplotypeFinderUnitTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/KBestHaplotypeFinderUnitTest.java
@@ -65,7 +65,6 @@ import java.util.*;
* User: rpoplin
* Date: 1/31/13
*/
-
public class KBestHaplotypeFinderUnitTest extends BaseTest {
@DataProvider(name = "BasicPathFindingData")
@@ -113,11 +112,11 @@ public class KBestHaplotypeFinderUnitTest extends BaseTest {
final int expectedNumOfPaths = nStartNodes * nBranchesPerBubble * nEndNodes;
Assert.assertEquals(paths.size(), expectedNumOfPaths, "Didn't find the expected number of paths");
- int lastScore = Integer.MAX_VALUE;
+ double lastScore = 0;
for ( final KBestHaplotype kbh : paths ) {
final Path path = kbh.path();
- Assert.assertTrue(path.getScore() <= lastScore, "Paths out of order. Path " + path + " has score above previous " + lastScore);
- lastScore = path.getScore();
+ Assert.assertTrue(kbh.score() <= lastScore, "Paths out of order. Path " + path + " has score " + path.getScore() + " above previous " + lastScore);
+ lastScore = kbh.score();
}
// get the best path, and make sure it's the same as our optimal path overall
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitterUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitterUnitTest.java
index 2f44129d8..eb0432769 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitterUnitTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/SharedVertexSequenceSplitterUnitTest.java
@@ -47,6 +47,7 @@
package org.broadinstitute.sting.gatk.walkers.haplotypecaller.graphs;
import org.broadinstitute.sting.BaseTest;
+import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.collections.Pair;
import org.testng.Assert;
@@ -226,28 +227,34 @@ public class SharedVertexSequenceSplitterUnitTest extends BaseTest {
}
final Set haplotypes = new HashSet<>();
- final List originalPaths = new KBestHaplotypeFinder((SeqGraph) graph.clone(),graph.getSources(),graph.getSinks());
+ final KBestHaplotypeFinder originalPaths = new KBestHaplotypeFinder((SeqGraph) graph.clone(),graph.getSources(),graph.getSinks());
for ( final KBestHaplotype path : originalPaths )
haplotypes.add(new String(path.bases()));
final SharedVertexSequenceSplitter splitter = new SharedVertexSequenceSplitter(graph, v);
splitter.split();
- if ( PRINT_GRAPHS ) graph.printGraph(new File(Utils.join("_", strings) + ".original.dot"), 0);
- if ( PRINT_GRAPHS ) splitter.splitGraph.printGraph(new File(Utils.join("_", strings) + ".split.dot"), 0);
+ if ( PRINT_GRAPHS ) graph.printGraph(new File(Utils.join("_", strings) + "_" + hasTop + "_" + hasBot + ".original.dot"), 0);
+ if ( PRINT_GRAPHS ) splitter.splitGraph.printGraph(new File(Utils.join("_", strings) + "_" + hasTop + "_" + hasBot + ".split.dot"), 0);
splitter.updateGraph(top, bot);
- if ( PRINT_GRAPHS ) graph.printGraph(new File(Utils.join("_", strings) + ".updated.dot"), 0);
+ if ( PRINT_GRAPHS ) graph.printGraph(new File(Utils.join("_", strings) + "_" + hasTop + "_" + hasBot + ".updated.dot"), 0);
- final List splitPaths = new KBestHaplotypeFinder(graph,graph.getSources(),graph.getSinks());
+ final KBestHaplotypeFinder splitPaths = new KBestHaplotypeFinder(graph,graph.getSources(),graph.getSinks());
for ( final KBestHaplotype path : splitPaths ) {
final String h = new String(path.bases());
Assert.assertTrue(haplotypes.contains(h), "Failed to find haplotype " + h);
}
- if ( splitPaths.size() == originalPaths.size() ) {
- for ( int i = 0; i < originalPaths.size(); i++ ) {
- Assert.assertTrue(splitPaths.get(i).path().equalScoreAndSequence(originalPaths.get(i).path()), "Paths not equal " + splitPaths.get(i) + " vs. original " + originalPaths.get(i));
- }
- }
+
+ final List sortedOriginalPaths = new ArrayList<>(originalPaths.size());
+ for (final KBestHaplotype kbh : originalPaths.unique())
+ sortedOriginalPaths.add(kbh.bases());
+ Collections.sort(sortedOriginalPaths, BaseUtils.BASES_COMPARATOR);
+ final List sortedSplitPaths = new ArrayList<>(splitPaths.size());
+ for (final KBestHaplotype kbh : splitPaths.unique())
+ sortedSplitPaths.add(kbh.bases());
+ Collections.sort(sortedSplitPaths, BaseUtils.BASES_COMPARATOR);
+
+ Assert.assertEquals(sortedSplitPaths,sortedOriginalPaths,Utils.join("_", strings) + "_" + hasTop + "_" + hasBot);
}
@DataProvider(name = "MeetsMinSequenceData")
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/indels/IndelRealignerIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/indels/IndelRealignerIntegrationTest.java
index 275ababda..185082c2d 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/indels/IndelRealignerIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/indels/IndelRealignerIntegrationTest.java
@@ -60,8 +60,8 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
private static final String knownIndels = validationDataLocation + "indelRealignerTest.pilot1.ceu.vcf";
private static final String baseCommandPrefix = "-T IndelRealigner -noPG -R " + b36KGReference + " -I " + mainTestBam + " -targetIntervals " + mainTestIntervals + " -compress 0 -L 20:49,500-55,500 ";
private static final String baseCommand = baseCommandPrefix + "-o %s ";
- private static final String base_md5 = "a102dd55451799e5f053c784b762087e";
- private static final String base_md5_with_SW_or_VCF = "06b8eefcbd785e929027feaa22bb060d";
+ private static final String base_md5 = "458588d68c8ea7e54443ea722604b265";
+ private static final String base_md5_with_SW_or_VCF = "d5ed91bd5b2023c69078a0fc00268d3c";
@Test
public void testDefaults() {
@@ -84,7 +84,7 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
WalkerTestSpec spec1 = new WalkerTestSpec(
baseCommand + "--consensusDeterminationModel KNOWNS_ONLY -known " + knownIndels,
1,
- Arrays.asList("1b24b0f2a20aed1adc726d1b296a3192"));
+ Arrays.asList("a1b9396f4d5b65f7ae6e0062daf363a3"));
executeTest("realigner known indels only from VCF", spec1);
}
@@ -101,7 +101,7 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
public void testLods() {
HashMap e = new HashMap();
e.put( "-LOD 60", base_md5 );
- e.put( "-LOD 1 --consensusDeterminationModel USE_SW", "4bf28d3c0337682d439257874377a681" );
+ e.put( "-LOD 1 --consensusDeterminationModel USE_SW", "dea9bd14323b33348d9cf28e256415f2" );
for ( Map.Entry entry : e.entrySet() ) {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
@@ -117,7 +117,7 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
WalkerTestSpec spec = new WalkerTestSpec(
"-T IndelRealigner -noPG -R " + b36KGReference + " -I " + validationDataLocation + "NA12878.chrom1.SLX.SRP000032.2009_06.bam -L 1:10,000,000-11,000,000 -targetIntervals " + validationDataLocation + "indelRealignerTest.NA12878.chrom1.intervals -compress 0 -o %s",
1,
- Arrays.asList("f4f6c3b2a2be0306a0ecd3def334bafe"));
+ Arrays.asList("b91c0bf803247f703dc1cb6ccdc4f18f"));
executeTest("realigner long run", spec);
}
@@ -126,7 +126,7 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
WalkerTestSpec spec = new WalkerTestSpec(
baseCommand + "--noOriginalAlignmentTags --consensusDeterminationModel USE_SW",
1,
- Arrays.asList("71fb521f8febfe2dc683fc636e28ae7d"));
+ Arrays.asList("041e2254f271261fb46dc3878cf638f6"));
executeTest("realigner no output tags", spec);
}
@@ -148,7 +148,7 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
@Test
public void testMaxReadsInMemory() {
HashMap e = new HashMap();
- e.put( "--maxReadsInMemory 10000", base_md5 );
+ e.put( "--maxReadsInMemory 10000", "0108cd5950f1a4eb90209c3dca8f9e11" );
e.put( "--maxReadsInMemory 40000", base_md5 );
for ( Map.Entry entry : e.entrySet() ) {
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFsIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFsIntegrationTest.java
index 03d136290..1a4d0c1f4 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFsIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/variantutils/CombineGVCFsIntegrationTest.java
@@ -165,4 +165,12 @@ public class CombineGVCFsIntegrationTest extends WalkerTest {
spec.disableShadowBCF();
executeTest("testMD5s", spec);
}
+
+ @Test
+ public void testBasepairResolution() throws Exception {
+ final String cmd = baseTestString(" -L 1:69485-69791 --convertToBasePairResolution");
+ final WalkerTestSpec spec = new WalkerTestSpec(cmd, 1, Arrays.asList("a068fb2c35cdd14df1e8f1f92f4114b4"));
+ spec.disableShadowBCF();
+ executeTest("testBasepairResolution", spec);
+ }
}
diff --git a/protected/pom.xml b/protected/pom.xml
index 8a9646438..434ea568e 100644
--- a/protected/pom.xml
+++ b/protected/pom.xml
@@ -5,7 +5,7 @@
org.broadinstitute.sting
sting-root
- 3.0
+ 3.1
../public/sting-root
diff --git a/public/VectorPairHMM/README.md b/public/VectorPairHMM/README.md
new file mode 100644
index 000000000..85cc0a04a
--- /dev/null
+++ b/public/VectorPairHMM/README.md
@@ -0,0 +1,71 @@
+Implementation overview:
+Created a new Java class called VectorLoglessPairHMM which extends LoglessPairHMM and
+overrides functions from both LoglessPairHMM and PairHMM.
+1. Constructor: Call base class constructors. Then, load the native library located in this
+directory and call an init function (with suffix 'jniInitializeClassFieldsAndMachineMask') in the
+library to determine fields ids for the members of classes JNIReadDataHolder and
+JNIHaplotypeDataHolders. The native code stores the field ids (struct offsets) for the classes and
+re-uses them for subsequent computations. Optionally, the user can disable the vector
+implementation, by using the 'mask' argument (see comments for a more detailed explanation).
+2. When the library is loaded, it invokes the constructor of the class LoadTimeInitializer (because
+a global variable g_load_time_initializer is declared in the library). This constructor
+(LoadTimeInitializer.cc) can be used to perform various initializations. Currently, it initializes
+two global function pointers to point to the function implementation that is supported on the
+machine (AVX/SSE/un-vectorized) on which the program is being run. The two pointers are for float
+and double respectively. The global function pointers are declared in utils.cc and are assigned in
+the function initialize_function_pointers() defined in utils.cc and invoked from the constructor of
+LoadTimeInitializer.
+Other initializations in LoadTimeInitializer:
+* ConvertChar::init - sets some masks for the vector implementation
+* FTZ for performance
+* stat counters = 0
+* debug structs (which are never used in non-debug mode)
+This initialization is done only once for the whole program.
+3. initialize(): To initialize the region for PairHMM. Pass haplotype bases to native code through
+the JNIHaplotypeDataHolder class. Since the haplotype list is common across multiple samples in
+computeReadLikelihoods(), we can pass the haplotype bases to the native code once and re-use across
+multiple samples.
+4. computeLikelihoods(): Copies array references for readBases/quals etc to array of
+JNIReadDataHolder objects. Invokes the JNI function to perform the computation and updates the
+likelihoodMap.
+The JNI function copies the byte array references into an array of testcase structs and invokes the
+compute_full_prob function through the function pointers initialized earlier.
+The primary native function called is
+Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniComputeLikelihoods. It uses
+standard JNI calls to get and return data from/to the Java class VectorLoglessPairHMM. The last
+argument to the function is the maximum number of OpenMP threads to use while computing PairHMM in
+C++. This option is set when the native function call is made from JNILoglessPairHMM
+computeLikelihoods - currently it is set to 12 (no logical reason).
+Note: OpenMP has been disabled for now - insufficient #testcases per call to computeLikelihoods() to
+justify multi-threading.
+5. finalizeRegion(): Releases the haplotype arrays initialized in step 3 - should be called at the
+end of every region (line 351 in PairHMMLikelihoodCalculationEngine).
+
+Note: Debug code has been moved to a separate class DebugJNILoglessPairHMM.java.
+
+Compiling:
+Make sure you have icc (Intel C compiler) available. Currently, gcc does not seem to support all AVX
+intrinsics.
+This native library is called libVectorLoglessPairHMM.so
+Using Maven:
+Type 'mvn install' in this directory - this will build the library (by invoking 'make') and copy the
+native library to the directory
+${sting-utils.basedir}/src/main/resources/org/broadinstitute/sting/utils/pairhmm
+The GATK maven build process (when run) will bundle the library into the StingUtils jar file from
+the copied directory.
+Simple build:
+cd src/main/c++
+make
+
+Running:
+The default implementation of PairHMM is now VECTOR_LOGLESS_CACHING in HaplotypeCaller.java. To use
+the Java version, use the command line argument "--pair_hmm_implementation LOGLESS_CACHING". (see
+run.sh in src/main/c++).
+The native library is bundled with the StingUtils jar file. When HaplotypeCaller is invoked, then
+the library is unpacked from the jar file, copied to the /tmp directory (with a unique id) and
+loaded by the Java class VectorLoglessPairHMM in the constructor (if it has not been loaded
+already).
+The default library can be overridden by using the -Djava.library.path argument (see
+src/main/c++/run.sh for an example) for the JVM to pass the path to the library. If the library
+libVectorLoglessPairHMM.so can be found in java.library.path, then it is loaded and the 'packed'
+library is not used.
diff --git a/public/VectorPairHMM/pom.xml b/public/VectorPairHMM/pom.xml
new file mode 100644
index 000000000..41bb73211
--- /dev/null
+++ b/public/VectorPairHMM/pom.xml
@@ -0,0 +1,119 @@
+
+
+ 4.0.0
+
+
+ org.broadinstitute.sting
+ sting-root
+ 2.8-SNAPSHOT
+ ../../public/sting-root
+
+
+ VectorPairHMM
+ pom
+ Vectorized PairHMM native libraries
+
+ Builds a GNU/Linux x86_64 library of VectorPairHMM using icc (Intel C++ compiler). During install, copies it into sting-utils. Neither tested nor expected to work on any other platform.
+
+
+ UTF-8
+ ${sourceEncoding}
+ ${sourceEncoding}
+ ${project.basedir}/../..
+ ${sting.basedir}/public/sting-utils
+
+ ${sting-utils.basedir}/src/main/resources/org/broadinstitute/sting/utils/pairhmm
+
+
+
+
+
+
+
+ org.apache.maven.plugins
+ maven-enforcer-plugin
+
+
+
+ display-info
+
+ validate
+
+
+
+
+
+
+ org.codehaus.mojo
+ exec-maven-plugin
+
+
+
+ exec
+
+ compile
+
+ make
+ src/main/c++
+
+ ${java.home}
+ ${project.build.directory}
+
+
+
+
+
+
+
+
+ org.apache.maven.plugins
+ maven-install-plugin
+
+ true
+
+
+
+
+
+ org.apache.maven.plugins
+ maven-resources-plugin
+
+
+ default-install
+
+ copy-resources
+
+ install
+
+ ${pairhmm.resources.directory}
+
+
+ ${project.build.directory}
+
+ **/*
+
+
+
+
+
+
+
+
+
+
+ com.google.code.sortpom
+ maven-sortpom-plugin
+
+ false
+ custom_1
+ \n
+ ${sourceEncoding}
+ true
+ scope
+ 4
+ false
+
+
+
+
+
diff --git a/public/VectorPairHMM/src/main/c++/.gitignore b/public/VectorPairHMM/src/main/c++/.gitignore
new file mode 100644
index 000000000..d791ffd80
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/.gitignore
@@ -0,0 +1,16 @@
+.svn
+*.o
+*.so
+tests
+.deps
+hmm_Mohammad
+pairhmm-template-main
+*.swp
+*.class
+checker
+reformat
+subdir_checkout.sh
+avx/
+sse/
+triplicate.sh
+
diff --git a/public/VectorPairHMM/src/main/c++/LoadTimeInitializer.cc b/public/VectorPairHMM/src/main/c++/LoadTimeInitializer.cc
new file mode 100644
index 000000000..0e3026f65
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/LoadTimeInitializer.cc
@@ -0,0 +1,206 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#include "LoadTimeInitializer.h"
+#include "utils.h"
+using namespace std;
+char* LoadTimeInitializerStatsNames[] =
+{
+ "num_regions",
+ "num_reads",
+ "num_haplotypes",
+ "num_testcases",
+ "num_double_invocations",
+ "haplotype_length",
+ "readlength",
+ "product_read_length_haplotype_length",
+ "dummy"
+};
+
+LoadTimeInitializer g_load_time_initializer;
+
+LoadTimeInitializer::LoadTimeInitializer() //will be called when library is loaded
+{
+ ConvertChar::init();
+#ifndef DISABLE_FTZ
+ //Very important to get good performance on Intel processors
+ //Function: enabling FTZ converts denormals to 0 in hardware
+ //Denormals cause microcode to insert uops into the core causing big slowdown
+ _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
+ cout << "FTZ enabled - may decrease accuracy if denormal numbers encountered\n";
+#else
+ cout << "FTZ is not set - may slow down performance if denormal numbers encountered\n";
+#endif
+ //Profiling: times for compute and transfer (either bytes copied or pointers copied)
+ m_compute_time = 0;
+ m_data_transfer_time = 0;
+ m_bytes_copied = 0;
+
+ //Initialize profiling counters
+ for(unsigned i=0;i::initializeStaticMembers();
+ Context::initializeStaticMembers();
+
+ cout.flush();
+}
+
+void LoadTimeInitializer::print_profiling()
+{
+ double mean = 0;
+ double variance = 0;
+ uint64_t denominator = 1;
+ cout << "Time spent in compute_testcases "< C++) "<::iterator mI = m_filename_to_fptr.find(filename);
+ ofstream* fptr = 0;
+ if(mI == m_filename_to_fptr.end())
+ {
+ m_filename_to_fptr[filename] = new ofstream();
+ fptr = m_filename_to_fptr[filename];
+ //File never seen before
+ if(m_written_files_set.find(filename) == m_written_files_set.end())
+ {
+ to_append = false;
+ m_written_files_set.insert(filename);
+ }
+ fptr->open(filename.c_str(), to_append ? ios::app : ios::out);
+ assert(fptr->is_open());
+ }
+ else
+ fptr = (*mI).second;
+ //ofstream fptr;
+ //fptr.open(filename.c_str(), to_append ? ofstream::app : ofstream::out);
+ (*fptr) << s;
+ if(add_newline)
+ (*fptr) << "\n";
+ //fptr.close();
+}
+void LoadTimeInitializer::debug_close()
+{
+ for(map::iterator mB = m_filename_to_fptr.begin(), mE = m_filename_to_fptr.end();
+ mB != mE;mB++)
+ {
+ (*mB).second->close();
+ delete (*mB).second;
+ }
+ m_filename_to_fptr.clear();
+}
+
+void LoadTimeInitializer::dump_sandbox(testcase& tc, unsigned tc_idx, unsigned numReads, unsigned numHaplotypes)
+{
+ unsigned haplotypeLength = tc.haplen;
+ unsigned readLength = tc.rslen;
+ ofstream& dumpFptr = m_sandbox_fptr;
+ for(unsigned k=0;k
+#include "template.h"
+
+enum LoadTimeInitializerStatsEnum
+{
+ NUM_REGIONS_IDX=0,
+ NUM_READS_IDX,
+ NUM_HAPLOTYPES_IDX,
+ NUM_TESTCASES_IDX,
+ NUM_DOUBLE_INVOCATIONS_IDX,
+ HAPLOTYPE_LENGTH_IDX,
+ READ_LENGTH_IDX,
+ PRODUCT_READ_LENGTH_HAPLOTYPE_LENGTH_IDX,
+ TOTAL_NUMBER_STATS
+};
+extern char* LoadTimeInitializerStatsNames[];
+
+class LoadTimeInitializer
+{
+ public:
+ LoadTimeInitializer(); //will be called when library is loaded
+ ~LoadTimeInitializer()
+ {
+ delete[] m_buffer;
+ }
+ void print_profiling();
+ void debug_dump(std::string filename, std::string s, bool to_append, bool add_newline=true);
+ void debug_close();
+
+ void dump_sandbox(testcase& tc, unsigned tc_idx, unsigned numReads, unsigned numHaplotypes);
+ void open_sandbox() { m_sandbox_fptr.open("sandbox.txt", std::ios::app); }
+ void close_sandbox() { m_sandbox_fptr.close(); }
+
+ jfieldID m_readBasesFID;
+ jfieldID m_readQualsFID;
+ jfieldID m_insertionGOPFID;
+ jfieldID m_deletionGOPFID;
+ jfieldID m_overallGCPFID;
+ jfieldID m_haplotypeBasesFID;
+ //profiling - update stats
+ void update_stat(LoadTimeInitializerStatsEnum stat_idx, uint64_t value);
+ //timing in nanoseconds
+ uint64_t m_compute_time;
+ uint64_t m_data_transfer_time;
+ //bytes copied
+ uint64_t m_bytes_copied;
+ unsigned get_buffer_size() { return m_buffer_size; }
+ char* get_buffer() { return (char*)m_buffer; }
+ private:
+ std::map m_filename_to_fptr;
+ std::set m_written_files_set;
+ std::ofstream m_sandbox_fptr;
+ //used to compute various stats
+ uint64_t m_sum_stats[TOTAL_NUMBER_STATS];
+ double m_sum_square_stats[TOTAL_NUMBER_STATS];
+ uint64_t m_min_stats[TOTAL_NUMBER_STATS];
+ uint64_t m_max_stats[TOTAL_NUMBER_STATS];
+ unsigned m_buffer_size;
+ uint64_t* m_buffer;
+};
+extern LoadTimeInitializer g_load_time_initializer;
+
+#define SIZE_PER_TESTCASE 6*10000
+#define SIZE_PER_BUFFER 10000
+
+#endif
diff --git a/public/VectorPairHMM/src/main/c++/Makefile b/public/VectorPairHMM/src/main/c++/Makefile
new file mode 100644
index 000000000..354bca0bb
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/Makefile
@@ -0,0 +1,114 @@
+#Copyright (c) 2012 The Broad Institute
+
+#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
+#NONINFRINGEMENT. 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.
+#
+
+
+#OMPCFLAGS=-fopenmp
+#OMPLFLAGS=-fopenmp #-openmp-link static
+
+#CFLAGS=-O2 -std=c++11 -W -Wall -march=corei7-avx -Wa,-q -pedantic $(OMPCFLAGS) -Wno-unknown-pragmas
+#CFLAGS=-O2 -W -Wall -march=corei7 -mfpmath=sse -msse4.2 -pedantic $(OMPCFLAGS) -Wno-unknown-pragmas
+
+JRE_HOME?=/opt/jdk1.7.0_25/jre
+JNI_COMPILATION_FLAGS=-D_REENTRANT -fPIC -I${JRE_HOME}/../include -I${JRE_HOME}/../include/linux
+
+COMMON_COMPILATION_FLAGS=$(JNI_COMPILATION_FLAGS) -O3 -W -Wall -pedantic $(OMPCFLAGS) -Wno-unknown-pragmas
+CC=icc
+CXX=icc
+
+LDFLAGS=-lm -lrt $(OMPLDFLAGS)
+ifdef DISABLE_FTZ
+ COMMON_COMPILATION_FLAGS+=-DDISABLE_FTZ -no-ftz
+endif
+
+PAPI_DIR=/home/karthikg/softwares/papi-5.3.0
+ifdef USE_PAPI
+ ifeq ($(USE_PAPI),1)
+ COMMON_COMPILATION_FLAGS+=-I$(PAPI_DIR)/include -DUSE_PAPI
+ LDFLAGS+=-L$(PAPI_DIR)/lib -lpapi
+ endif
+endif
+
+ifdef DISABLE_FTZ
+ COMMON_COMPILATION_FLAGS+=-DDISABLE_FTZ -no-ftz
+endif
+
+BIN=libVectorLoglessPairHMM.so pairhmm-template-main checker
+#BIN=checker
+
+DEPDIR=.deps
+DF=$(DEPDIR)/$(*).d
+
+#Common across libJNI and sandbox
+COMMON_SOURCES=utils.cc avx_function_instantiations.cc baseline.cc sse_function_instantiations.cc LoadTimeInitializer.cc
+#Part of libJNI
+LIBSOURCES=org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc Sandbox.cc $(COMMON_SOURCES)
+SOURCES=$(LIBSOURCES) pairhmm-template-main.cc pairhmm-1-base.cc
+LIBOBJECTS=$(LIBSOURCES:.cc=.o)
+COMMON_OBJECTS=$(COMMON_SOURCES:.cc=.o)
+
+
+#No vectorization for these files
+NO_VECTOR_SOURCES=org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc pairhmm-template-main.cc pairhmm-1-base.cc utils.cc baseline.cc LoadTimeInitializer.cc Sandbox.cc
+#Use -xAVX for these files
+AVX_SOURCES=avx_function_instantiations.cc
+#Use -xSSE4.2 for these files
+SSE_SOURCES=sse_function_instantiations.cc
+
+NO_VECTOR_OBJECTS=$(NO_VECTOR_SOURCES:.cc=.o)
+AVX_OBJECTS=$(AVX_SOURCES:.cc=.o)
+SSE_OBJECTS=$(SSE_SOURCES:.cc=.o)
+$(NO_VECTOR_OBJECTS): CXXFLAGS=$(COMMON_COMPILATION_FLAGS)
+$(AVX_OBJECTS): CXXFLAGS=$(COMMON_COMPILATION_FLAGS) -xAVX
+$(SSE_OBJECTS): CXXFLAGS=$(COMMON_COMPILATION_FLAGS) -xSSE4.2
+OBJECTS=$(NO_VECTOR_OBJECTS) $(AVX_OBJECTS) $(SSE_OBJECTS)
+
+all: $(BIN) Sandbox.class copied_lib
+
+-include $(addprefix $(DEPDIR)/,$(SOURCES:.cc=.d))
+
+checker: pairhmm-1-base.o $(COMMON_OBJECTS)
+ $(CXX) $(OMPLFLAGS) -o $@ $^ $(LDFLAGS)
+
+pairhmm-template-main: pairhmm-template-main.o $(COMMON_OBJECTS)
+ $(CXX) $(OMPLFLAGS) -o $@ $^ $(LDFLAGS)
+
+libVectorLoglessPairHMM.so: $(LIBOBJECTS)
+ $(CXX) $(OMPLFLAGS) -shared -static-intel -o $@ $(LIBOBJECTS) ${LDFLAGS}
+
+
+$(OBJECTS): %.o: %.cc
+ @mkdir -p $(DEPDIR)
+ $(CXX) -c -MMD -MF $(DF) $(CXXFLAGS) $(OUTPUT_OPTION) $<
+
+Sandbox.class: Sandbox.java
+ javac Sandbox.java
+
+copied_lib: libVectorLoglessPairHMM.so
+ifdef OUTPUT_DIR
+ mkdir -p $(OUTPUT_DIR)
+ rsync -a libVectorLoglessPairHMM.so $(OUTPUT_DIR)/
+endif
+
+clean:
+ rm -rf $(BIN) *.o $(DEPDIR) *.class
diff --git a/public/VectorPairHMM/src/main/c++/Sandbox.cc b/public/VectorPairHMM/src/main/c++/Sandbox.cc
new file mode 100644
index 000000000..985b19ae9
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/Sandbox.cc
@@ -0,0 +1,106 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#include "Sandbox.h"
+#include "org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.h"
+#include "utils.h"
+#include "jni_common.h"
+/*
+ * Class: Sandbox
+ * Method: jniGetMachineType
+ * Signature: ()J
+ */
+JNIEXPORT jlong JNICALL Java_Sandbox_jniGetMachineType
+ (JNIEnv * env, jobject thisObj)
+{
+ return 0;
+}
+
+/*
+ * Class: Sandbox
+ * Method: jniInitializeClassFieldsAndMachineMask
+ * Signature: (Ljava/lang/Class;Ljava/lang/Class;J)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniInitializeClassFieldsAndMachineMask
+ (JNIEnv* env, jobject thisObject, jclass readDataHolderClass, jclass haplotypeDataHolderClass, jlong mask)
+{
+ Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniInitializeClassFieldsAndMachineMask(env, thisObject, readDataHolderClass,
+ haplotypeDataHolderClass, mask);
+}
+
+/*
+ * Class: Sandbox
+ * Method: jniInitializeHaplotypes
+ * Signature: (I[LSandbox/JNIHaplotypeDataHolderClass;)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniInitializeHaplotypes
+ (JNIEnv * env, jobject thisObject, jint numHaplotypes, jobjectArray haplotypeDataArray)
+{
+ Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniInitializeHaplotypes(env, thisObject, numHaplotypes, haplotypeDataArray);
+}
+
+/*
+ * Class: Sandbox
+ * Method: jniFinalizeRegion
+ * Signature: ()V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniFinalizeRegion
+ (JNIEnv * env, jobject thisObject)
+{
+ Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniFinalizeRegion(env, thisObject);
+}
+
+
+/*
+ * Class: Sandbox
+ * Method: jniComputeLikelihoods
+ * Signature: (II[LSandbox/JNIReadDataHolderClass;[LSandbox/JNIHaplotypeDataHolderClass;[DI)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniComputeLikelihoods
+ (JNIEnv* env, jobject thisObject, jint numReads, jint numHaplotypes,
+ jobjectArray readDataArray, jobjectArray haplotypeDataArray, jdoubleArray likelihoodArray, jint maxNumThreadsToUse)
+{
+ Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniComputeLikelihoods(env, thisObject,
+ numReads, numHaplotypes, readDataArray, haplotypeDataArray, likelihoodArray, maxNumThreadsToUse);
+}
+/*
+ * Class: Sandbox
+ * Method: jniClose
+ * Signature: ()V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniClose
+ (JNIEnv* env, jobject thisObject)
+{ Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniClose(env, thisObject); }
+
+JNIEXPORT void JNICALL Java_Sandbox_doEverythingNative
+ (JNIEnv* env, jobject thisObject, jstring fileNameString)
+{
+ const char* fileName = env->GetStringUTFChars(fileNameString, 0);
+ char local_array[800];
+ strncpy(local_array, fileName, 200);
+ env->ReleaseStringUTFChars(fileNameString, fileName);
+ do_compute(local_array, true, 10000, false);
+}
+
diff --git a/public/VectorPairHMM/src/main/c++/Sandbox.h b/public/VectorPairHMM/src/main/c++/Sandbox.h
new file mode 100644
index 000000000..486a1c095
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/Sandbox.h
@@ -0,0 +1,96 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+/* DO NOT EDIT THIS FILE - it is machine generated */
+#include
+/* Header for class Sandbox */
+
+#ifndef _Included_Sandbox
+#define _Included_Sandbox
+#ifdef __cplusplus
+extern "C" {
+#endif
+#undef Sandbox_enableAll
+#define Sandbox_enableAll -1LL
+/*
+ * Class: Sandbox
+ * Method: jniGetMachineType
+ * Signature: ()J
+ */
+JNIEXPORT jlong JNICALL Java_Sandbox_jniGetMachineType
+ (JNIEnv *, jobject);
+
+/*
+ * Class: Sandbox
+ * Method: jniInitializeClassFieldsAndMachineMask
+ * Signature: (Ljava/lang/Class;Ljava/lang/Class;J)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniInitializeClassFieldsAndMachineMask
+ (JNIEnv *, jobject, jclass, jclass, jlong);
+
+/*
+ * Class: Sandbox
+ * Method: jniInitializeHaplotypes
+ * Signature: (I[LSandbox/JNIHaplotypeDataHolderClass;)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniInitializeHaplotypes
+ (JNIEnv *, jobject, jint, jobjectArray);
+
+/*
+ * Class: Sandbox
+ * Method: jniFinalizeRegion
+ * Signature: ()V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniFinalizeRegion
+ (JNIEnv *, jobject);
+
+/*
+ * Class: Sandbox
+ * Method: jniComputeLikelihoods
+ * Signature: (II[LSandbox/JNIReadDataHolderClass;[LSandbox/JNIHaplotypeDataHolderClass;[DI)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniComputeLikelihoods
+ (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jdoubleArray, jint);
+
+/*
+ * Class: Sandbox
+ * Method: jniClose
+ * Signature: ()V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_jniClose
+ (JNIEnv *, jobject);
+
+/*
+ * Class: Sandbox
+ * Method: doEverythingNative
+ * Signature: ([B)V
+ */
+JNIEXPORT void JNICALL Java_Sandbox_doEverythingNative
+ (JNIEnv *, jobject, jstring);
+
+#ifdef __cplusplus
+}
+#endif
+#endif
diff --git a/public/VectorPairHMM/src/main/c++/Sandbox.java b/public/VectorPairHMM/src/main/c++/Sandbox.java
new file mode 100644
index 000000000..698c1b48c
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/Sandbox.java
@@ -0,0 +1,306 @@
+/*
+* Copyright (c) 2012 The Broad Institute
+*
+* 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
+* NONINFRINGEMENT. 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 org.broadinstitute.sting.utils.vectorpairhmm;
+
+import java.util.List;
+import java.util.LinkedList;
+import java.util.Map;
+import java.util.HashMap;
+import java.io.File;
+import java.util.Scanner;
+import java.io.IOException;
+import java.io.FileNotFoundException;
+import java.io.InputStreamReader;
+
+public class Sandbox {
+
+ private long setupTime = 0;
+ private long computeTime = 0;
+ //Used to copy references to byteArrays to JNI from reads
+ protected class JNIReadDataHolderClass {
+ public byte[] readBases = null;
+ public byte[] readQuals = null;
+ public byte[] insertionGOP = null;
+ public byte[] deletionGOP = null;
+ public byte[] overallGCP = null;
+ }
+
+ //Used to copy references to byteArrays to JNI from haplotypes
+ protected class JNIHaplotypeDataHolderClass {
+ public byte[] haplotypeBases = null;
+ }
+
+ /**
+ * Return 64-bit mask representing machine capabilities
+ * Bit 0 is LSB, bit 63 MSB
+ * Bit 0 represents sse4.2 availability
+ * Bit 1 represents AVX availability
+ */
+ public native long jniGetMachineType();
+ public static final long enableAll = 0xFFFFFFFFFFFFFFFFl;
+
+
+ /**
+ * Function to initialize the fields of JNIReadDataHolderClass and JNIHaplotypeDataHolderClass from JVM.
+ * C++ codegets FieldIDs for these classes once and re-uses these IDs for the remainder of the program. Field IDs do not
+ * change per JVM session
+ * @param readDataHolderClass class type of JNIReadDataHolderClass
+ * @param haplotypeDataHolderClass class type of JNIHaplotypeDataHolderClass
+ * @param mask mask is a 64 bit integer identical to the one received from jniGetMachineType(). Users can disable usage of some hardware features by zeroing some bits in the mask
+ * */
+ private native void jniInitializeClassFieldsAndMachineMask(Class readDataHolderClass, Class haplotypeDataHolderClass, long mask);
+
+ private static Boolean isVectorLoglessPairHMMLibraryLoaded = false;
+ //The constructor is called only once inside PairHMMLikelihoodCalculationEngine
+ public Sandbox() {
+ synchronized(isVectorLoglessPairHMMLibraryLoaded) {
+ //Load the library and initialize the FieldIDs
+ if(!isVectorLoglessPairHMMLibraryLoaded) {
+ System.loadLibrary("VectorLoglessPairHMM");
+ isVectorLoglessPairHMMLibraryLoaded = true;
+ jniInitializeClassFieldsAndMachineMask(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class, enableAll); //need to do this only once
+ }
+ }
+ }
+
+ private native void jniInitializeHaplotypes(final int numHaplotypes, JNIHaplotypeDataHolderClass[] haplotypeDataArray);
+ private JNIHaplotypeDataHolderClass[] mHaplotypeDataArray = null;
+
+ //Used to transfer data to JNI
+ //Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region
+ public void initialize(final List haplotypes) {
+ int numHaplotypes = haplotypes.size();
+ mHaplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes];
+ int idx = 0;
+ for(final JNIHaplotypeDataHolderClass currHaplotype : haplotypes)
+ {
+ mHaplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass();
+ mHaplotypeDataArray[idx].haplotypeBases = currHaplotype.haplotypeBases;
+ ++idx;
+ }
+ jniInitializeHaplotypes(numHaplotypes, mHaplotypeDataArray);
+ }
+ /**
+ * Tell JNI to release arrays - really important if native code is directly accessing Java memory, if not
+ * accessing Java memory directly, still important to release memory from C++
+ */
+ private native void jniFinalizeRegion();
+
+
+ public void finalizeRegion()
+ {
+ jniFinalizeRegion();
+ }
+
+ /**
+ * Real compute kernel
+ */
+ private native void jniComputeLikelihoods(int numReads, int numHaplotypes, JNIReadDataHolderClass[] readDataArray,
+ JNIHaplotypeDataHolderClass[] haplotypeDataArray, double[] likelihoodArray, int maxNumThreadsToUse);
+
+ public void computeLikelihoods(final List reads, final List haplotypes) {
+ //System.out.println("Region : "+reads.size()+" x "+haplotypes.size());
+ long startTime = System.nanoTime();
+ int readListSize = reads.size();
+ int numHaplotypes = haplotypes.size();
+ int numTestcases = readListSize*numHaplotypes;
+ JNIReadDataHolderClass[] readDataArray = new JNIReadDataHolderClass[readListSize];
+ int idx = 0;
+ for(JNIReadDataHolderClass read : reads)
+ {
+ readDataArray[idx] = new JNIReadDataHolderClass();
+ readDataArray[idx].readBases = read.readBases;
+ readDataArray[idx].readQuals = read.readQuals;
+ readDataArray[idx].insertionGOP = read.insertionGOP;
+ readDataArray[idx].deletionGOP = read.deletionGOP;
+ readDataArray[idx].overallGCP = read.overallGCP;
+ ++idx;
+ }
+
+ double[] mLikelihoodArray = new double[readListSize*numHaplotypes]; //to store results
+ setupTime += (System.nanoTime() - startTime);
+ //for(reads)
+ // for(haplotypes)
+ // compute_full_prob()
+ jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, mHaplotypeDataArray, mLikelihoodArray, 12);
+
+ computeTime += (System.nanoTime() - startTime);
+ }
+
+ /**
+ * Print final profiling information from native code
+ */
+ public native void jniClose();
+ public void close()
+ {
+ System.out.println("Time spent in setup for JNI call : "+(setupTime*1e-9)+" compute time : "+(computeTime*1e-9));
+ jniClose();
+ }
+
+ public void parseSandboxFile(String filename)
+ {
+ File file = new File(filename);
+ Scanner input = null;
+ try
+ {
+ input = new Scanner(file);
+ }
+ catch(FileNotFoundException e)
+ {
+ System.err.println("File "+filename+" cannot be found/read");
+ return;
+ }
+ int idx = 0;
+ int numReads = 0;
+ int numHaplotypes = 0;
+ int readIdx = 0, testCaseIdx = 0, haplotypeIdx = 0;
+ LinkedList haplotypeList = new LinkedList();
+ LinkedList readList = new LinkedList();
+
+ byte[][] byteArray = new byte[6][];
+ boolean firstLine = true;
+ String[] currTokens = new String[8];
+ while(input.hasNextLine())
+ {
+ String line = input.nextLine();
+ Scanner lineScanner = new Scanner(line);
+ idx = 0;
+ while(lineScanner.hasNext())
+ currTokens[idx++] = lineScanner.next();
+ if(idx == 0)
+ break;
+ assert(idx >= 6);
+ //start of new region
+ if(idx == 8)
+ {
+ if(!firstLine)
+ {
+ initialize(haplotypeList);
+ computeLikelihoods(readList, haplotypeList);
+ finalizeRegion();
+ }
+ try
+ {
+ numReads = Integer.parseInt(currTokens[6]);
+ }
+ catch(NumberFormatException e)
+ {
+ numReads = 1;
+ }
+ try
+ {
+ numHaplotypes = Integer.parseInt(currTokens[7]);
+ }
+ catch(NumberFormatException e)
+ {
+ numHaplotypes = 1;
+ }
+ haplotypeIdx = readIdx = testCaseIdx = 0;
+ readList.clear();
+ haplotypeList.clear();
+ }
+ if(haplotypeIdx < numHaplotypes)
+ {
+ JNIHaplotypeDataHolderClass X = new JNIHaplotypeDataHolderClass();
+ X.haplotypeBases = currTokens[0].getBytes();
+ haplotypeList.add(X);
+ }
+ if(testCaseIdx%numHaplotypes == 0)
+ {
+ JNIReadDataHolderClass X = new JNIReadDataHolderClass();
+ X.readBases = currTokens[1].getBytes();
+ for(int i=2;i<6;++i)
+ {
+ byteArray[i] = currTokens[i].getBytes();
+ for(int j=0;j 0 && readList.size() > 0)
+ {
+ initialize(haplotypeList);
+ computeLikelihoods(readList, haplotypeList);
+ finalizeRegion();
+ }
+
+ close();
+ input.close();
+ }
+
+ private native void doEverythingNative(String filename);
+
+ public static void main(String[] args)
+ {
+ if(args.length <= 0)
+ {
+ System.err.println("Needs 1 argument - ");
+ System.exit(-1);
+ }
+ //// Get runtime
+ //java.lang.Runtime rt = java.lang.Runtime.getRuntime();
+ //// Start a new process: UNIX command ls
+ //String cmd = "/home/karthikg/broad/gsa-unstable/public/c++/VectorPairHMM/checker "+args[0];
+ //try
+ //{
+ //System.out.println(cmd);
+ //java.lang.Process p = rt.exec(cmd);
+ //try
+ //{
+ //p.waitFor();
+ //java.io.InputStream is = p.getInputStream();
+ //java.io.BufferedReader reader = new java.io.BufferedReader(new InputStreamReader(is));
+ //// And print each line
+ //String s = null;
+ //while ((s = reader.readLine()) != null) {
+ //System.out.println(s);
+ //}
+ //is.close();
+ //}
+ //catch(InterruptedException e)
+ //{
+ //System.err.println(e);
+ //}
+ //}
+ //catch(IOException e)
+ //{
+ //System.err.println(e);
+ //}
+ Sandbox t = new Sandbox();
+ //t.doEverythingNative(args[0]);
+ t.parseSandboxFile(args[0]);
+ }
+}
diff --git a/public/VectorPairHMM/src/main/c++/Sandbox_JNIHaplotypeDataHolderClass.h b/public/VectorPairHMM/src/main/c++/Sandbox_JNIHaplotypeDataHolderClass.h
new file mode 100644
index 000000000..7f78f0178
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/Sandbox_JNIHaplotypeDataHolderClass.h
@@ -0,0 +1,13 @@
+/* DO NOT EDIT THIS FILE - it is machine generated */
+#include
+/* Header for class Sandbox_JNIHaplotypeDataHolderClass */
+
+#ifndef _Included_Sandbox_JNIHaplotypeDataHolderClass
+#define _Included_Sandbox_JNIHaplotypeDataHolderClass
+#ifdef __cplusplus
+extern "C" {
+#endif
+#ifdef __cplusplus
+}
+#endif
+#endif
diff --git a/public/VectorPairHMM/src/main/c++/Sandbox_JNIReadDataHolderClass.h b/public/VectorPairHMM/src/main/c++/Sandbox_JNIReadDataHolderClass.h
new file mode 100644
index 000000000..a9312ff3b
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/Sandbox_JNIReadDataHolderClass.h
@@ -0,0 +1,13 @@
+/* DO NOT EDIT THIS FILE - it is machine generated */
+#include
+/* Header for class Sandbox_JNIReadDataHolderClass */
+
+#ifndef _Included_Sandbox_JNIReadDataHolderClass
+#define _Included_Sandbox_JNIReadDataHolderClass
+#ifdef __cplusplus
+extern "C" {
+#endif
+#ifdef __cplusplus
+}
+#endif
+#endif
diff --git a/public/VectorPairHMM/src/main/c++/avx_function_instantiations.cc b/public/VectorPairHMM/src/main/c++/avx_function_instantiations.cc
new file mode 100644
index 000000000..6d90d5070
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/avx_function_instantiations.cc
@@ -0,0 +1,44 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#include "template.h"
+
+#undef SIMD_ENGINE
+#undef SIMD_ENGINE_SSE
+
+#define SIMD_ENGINE avx
+#define SIMD_ENGINE_AVX
+
+#include "define-float.h"
+#include "shift_template.c"
+#include "pairhmm-template-kernel.cc"
+
+#include "define-double.h"
+#include "shift_template.c"
+#include "pairhmm-template-kernel.cc"
+
+template double compute_full_prob_avxd(testcase* tc, double* nextlog);
+template float compute_full_prob_avxs(testcase* tc, float* nextlog);
+
diff --git a/public/VectorPairHMM/src/main/c++/baseline.cc b/public/VectorPairHMM/src/main/c++/baseline.cc
new file mode 100644
index 000000000..a2dd8d329
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/baseline.cc
@@ -0,0 +1,167 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#include "headers.h"
+#include "template.h"
+#include "utils.h"
+#include "LoadTimeInitializer.h"
+using namespace std;
+
+template
+NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log)
+{
+ int r, c;
+ int ROWS = tc->rslen + 1;
+ int COLS = tc->haplen + 1;
+
+ Context ctx;
+ //#define USE_STACK_ALLOCATION 1
+#ifdef USE_STACK_ALLOCATION
+ NUMBER M[ROWS][COLS];
+ NUMBER X[ROWS][COLS];
+ NUMBER Y[ROWS][COLS];
+ NUMBER p[ROWS][6];
+#else
+ //allocate on heap in way that simulates a 2D array. Having a 2D array instead of
+ //a straightforward array of pointers ensures that all data lies 'close' in memory, increasing
+ //the chance of being stored together in the cache. Also, prefetchers can learn memory access
+ //patterns for 2D arrays, not possible for array of pointers
+ //bool locally_allocated = false;
+ //NUMBER* common_buffer = 0;
+ NUMBER* common_buffer = new NUMBER[3*ROWS*COLS + ROWS*6];
+ //unsigned curr_size = sizeof(NUMBER)*(3*ROWS*COLS + ROWS*6);
+ //if(true)
+ //{
+ //common_buffer = new NUMBER[3*ROWS*COLS + ROWS*6];
+ //locally_allocated = true;
+ //}
+ //else
+ //common_buffer = (NUMBER*)(g_load_time_initializer.get_buffer());
+ //pointers to within the allocated buffer
+ NUMBER** common_pointer_buffer = new NUMBER*[4*ROWS];
+ NUMBER* ptr = common_buffer;
+ unsigned i = 0;
+ for(i=0;i<3*ROWS;++i, ptr+=COLS)
+ common_pointer_buffer[i] = ptr;
+ for(;i<4*ROWS;++i, ptr+=6)
+ common_pointer_buffer[i] = ptr;
+
+ NUMBER** M = common_pointer_buffer;
+ NUMBER** X = M + ROWS;
+ NUMBER** Y = X + ROWS;
+ NUMBER** p = Y + ROWS;
+#endif
+
+
+ p[0][MM] = ctx._(0.0);
+ p[0][GapM] = ctx._(0.0);
+ p[0][MX] = ctx._(0.0);
+ p[0][XX] = ctx._(0.0);
+ p[0][MY] = ctx._(0.0);
+ p[0][YY] = ctx._(0.0);
+
+ for (r = 1; r < ROWS; r++)
+ {
+ int _i = tc->i[r-1] & 127;
+ int _d = tc->d[r-1] & 127;
+ int _c = tc->c[r-1] & 127;
+ //p[r][MM] = ctx._(1.0) - ctx.ph2pr[(_i + _d) & 127];
+ SET_MATCH_TO_MATCH_PROB(p[r][MM], _i, _d);
+ p[r][GapM] = ctx._(1.0) - ctx.ph2pr[_c];
+ p[r][MX] = ctx.ph2pr[_i];
+ p[r][XX] = ctx.ph2pr[_c];
+ p[r][MY] = ctx.ph2pr[_d];
+ p[r][YY] = ctx.ph2pr[_c];
+ //p[r][MY] = (r == ROWS - 1) ? ctx._(1.0) : ctx.ph2pr[_d];
+ //p[r][YY] = (r == ROWS - 1) ? ctx._(1.0) : ctx.ph2pr[_c];
+ }
+ for (c = 0; c < COLS; c++)
+ {
+ M[0][c] = ctx._(0.0);
+ X[0][c] = ctx._(0.0);
+ Y[0][c] = ctx.INITIAL_CONSTANT / (tc->haplen);
+ }
+
+ for (r = 1; r < ROWS; r++)
+ {
+ M[r][0] = ctx._(0.0);
+ X[r][0] = X[r-1][0] * p[r][XX];
+ Y[r][0] = ctx._(0.0);
+ }
+
+ NUMBER result = ctx._(0.0);
+
+ for (r = 1; r < ROWS; r++)
+ for (c = 1; c < COLS; c++)
+ {
+ fexcept_t flagp;
+ char _rs = tc->rs[r-1];
+ char _hap = tc->hap[c-1];
+ int _q = tc->q[r-1] & 127;
+ NUMBER distm = ctx.ph2pr[_q];
+ if (_rs == _hap || _rs == 'N' || _hap == 'N')
+ distm = ctx._(1.0) - distm;
+ else
+ distm = distm/3;
+
+
+ //feclearexcept(FE_ALL_EXCEPT);
+ M[r][c] = distm * (M[r-1][c-1] * p[r][MM] + X[r-1][c-1] * p[r][GapM] + Y[r-1][c-1] * p[r][GapM]);
+ //STORE_FP_EXCEPTIONS(flagp, exceptions_array);
+
+ //feclearexcept(FE_ALL_EXCEPT);
+ X[r][c] = M[r-1][c] * p[r][MX] + X[r-1][c] * p[r][XX];
+ //STORE_FP_EXCEPTIONS(flagp, exceptions_array);
+
+ //feclearexcept(FE_ALL_EXCEPT);
+ Y[r][c] = M[r][c-1] * p[r][MY] + Y[r][c-1] * p[r][YY];
+ //STORE_FP_EXCEPTIONS(flagp, exceptions_array);
+
+ //CONVERT_AND_PRINT(M[r][c]);
+ //CONVERT_AND_PRINT(X[r][c]);
+ //CONVERT_AND_PRINT(Y[r][c]);
+
+ }
+ for (c = 0; c < COLS; c++)
+ {
+ result += M[ROWS-1][c] + X[ROWS-1][c];
+ }
+
+ if (before_last_log != NULL)
+ *before_last_log = result;
+
+#ifndef USE_STACK_ALLOCATION
+ delete[] common_pointer_buffer;
+ //if(locally_allocated)
+ delete[] common_buffer;
+#endif
+
+ return result;
+ //return ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT;
+}
+
+template double compute_full_prob(testcase* tc, double* nextbuf);
+template float compute_full_prob(testcase* tc, float* nextbuf);
+
diff --git a/public/VectorPairHMM/src/main/c++/define-double.h b/public/VectorPairHMM/src/main/c++/define-double.h
new file mode 100644
index 000000000..2067d369c
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/define-double.h
@@ -0,0 +1,205 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#include
+
+#ifdef PRECISION
+#undef PRECISION
+#undef MAIN_TYPE
+#undef MAIN_TYPE_SIZE
+#undef UNION_TYPE
+#undef IF_128
+#undef IF_MAIN_TYPE
+#undef SHIFT_CONST1
+#undef SHIFT_CONST2
+#undef SHIFT_CONST3
+#undef _128_TYPE
+#undef SIMD_TYPE
+#undef AVX_LENGTH
+#undef HAP_TYPE
+#undef MASK_TYPE
+#undef MASK_ALL_ONES
+
+#undef SET_VEC_ZERO(__vec)
+#undef VEC_OR(__v1, __v2)
+#undef VEC_ADD(__v1, __v2)
+#undef VEC_SUB(__v1, __v2)
+#undef VEC_MUL(__v1, __v2)
+#undef VEC_DIV(__v1, __v2)
+#undef VEC_BLEND(__v1, __v2, __mask)
+#undef VEC_BLENDV(__v1, __v2, __maskV)
+#undef VEC_CAST_256_128(__v1)
+#undef VEC_EXTRACT_128(__v1, __im)
+#undef VEC_EXTRACT_UNIT(__v1, __im)
+#undef VEC_SET1_VAL128(__val)
+#undef VEC_MOVE(__v1, __val)
+#undef VEC_CAST_128_256(__v1)
+#undef VEC_INSERT_VAL(__v1, __val, __pos)
+#undef VEC_CVT_128_256(__v1)
+#undef VEC_SET1_VAL(__val)
+#undef VEC_POPCVT_CHAR(__ch)
+#undef VEC_LDPOPCVT_CHAR(__addr)
+#undef VEC_CMP_EQ(__v1, __v2)
+#undef VEC_SET_LSE(__val)
+#undef SHIFT_HAP(__v1, __val)
+#undef MASK_VEC
+#undef VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst)
+#undef VEC_SHIFT_LEFT_1BIT(__vs)
+#undef MASK_ALL_ONES
+#undef COMPARE_VECS(__v1, __v2)
+#undef _256_INT_TYPE
+#undef BITMASK_VEC
+#endif
+
+#define PRECISION d
+#define MAIN_TYPE double
+#define MAIN_TYPE_SIZE 64
+#define UNION_TYPE mix_D
+#define IF_128 IF_128d
+#define IF_MAIN_TYPE IF_64
+#define SHIFT_CONST1 8
+#define SHIFT_CONST2 1
+#define SHIFT_CONST3 8
+#define _128_TYPE __m128d
+#define SIMD_TYPE __m256d
+#define _256_INT_TYPE __m256i
+#define AVX_LENGTH 4
+#define HAP_TYPE __m128i
+#define MASK_TYPE uint64_t
+#define MASK_ALL_ONES 0xFFFFFFFFFFFFFFFF
+#define MASK_VEC MaskVec_D
+
+#define SET_VEC_ZERO(__vec) \
+ __vec= _mm256_setzero_pd()
+
+#define VEC_OR(__v1, __v2) \
+ _mm256_or_pd(__v1, __v2)
+
+#define VEC_ADD(__v1, __v2) \
+ _mm256_add_pd(__v1, __v2)
+
+#define VEC_SUB(__v1, __v2) \
+ _mm256_sub_pd(__v1, __v2)
+
+#define VEC_MUL(__v1, __v2) \
+ _mm256_mul_pd(__v1, __v2)
+
+#define VEC_DIV(__v1, __v2) \
+ _mm256_div_pd(__v1, __v2)
+
+#define VEC_BLEND(__v1, __v2, __mask) \
+ _mm256_blend_pd(__v1, __v2, __mask)
+
+#define VEC_BLENDV(__v1, __v2, __maskV) \
+ _mm256_blendv_pd(__v1, __v2, __maskV)
+
+#define VEC_CAST_256_128(__v1) \
+ _mm256_castpd256_pd128 (__v1)
+
+#define VEC_EXTRACT_128(__v1, __im) \
+ _mm256_extractf128_pd (__v1, __im)
+
+#define VEC_EXTRACT_UNIT(__v1, __im) \
+ _mm_extract_epi64(__v1, __im)
+
+#define VEC_SET1_VAL128(__val) \
+ _mm_set1_pd(__val)
+
+#define VEC_MOVE(__v1, __val) \
+ _mm_move_sd(__v1, __val)
+
+#define VEC_CAST_128_256(__v1) \
+ _mm256_castpd128_pd256(__v1)
+
+#define VEC_INSERT_VAL(__v1, __val, __pos) \
+ _mm256_insertf128_pd(__v1, __val, __pos)
+
+#define VEC_CVT_128_256(__v1) \
+ _mm256_cvtepi32_pd(__v1)
+
+#define VEC_SET1_VAL(__val) \
+ _mm256_set1_pd(__val)
+
+#define VEC_POPCVT_CHAR(__ch) \
+ _mm256_cvtepi32_pd(_mm_set1_epi32(__ch))
+
+#define VEC_LDPOPCVT_CHAR(__addr) \
+ _mm256_cvtepi32_pd(_mm_load_si128((__m128i const *)__addr))
+
+#define VEC_CMP_EQ(__v1, __v2) \
+ _mm256_cmp_pd(__v1, __v2, _CMP_EQ_OQ)
+
+#define VEC_SET_LSE(__val) \
+ _mm256_set_pd(zero, zero, zero, __val);
+
+#define SHIFT_HAP(__v1, __val) \
+ __v1 = _mm_insert_epi32(_mm_slli_si128(__v1, 4), __val.i, 0)
+
+#define VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst) \
+ __vdst = _mm256_castpd128_pd256(__vsLow) ; \
+__vdst = _mm256_insertf128_pd(__vdst, __vsHigh, 1) ;
+
+#define VEC_SHIFT_LEFT_1BIT(__vs) \
+ __vs = _mm_slli_epi64(__vs, 1)
+
+
+#define COMPARE_VECS(__v1, __v2, __first, __last) { \
+ double* ptr1 = (double*) (&__v1) ; \
+ double* ptr2 = (double*) (&__v2) ; \
+ for (int ei=__first; ei <= __last; ++ei) { \
+ if (ptr1[ei] != ptr2[ei]) { \
+ std::cout << "Double Mismatch at " << ei << ": " \
+ << ptr1[ei] << " vs. " << ptr2[ei] << std::endl ; \
+ exit(0) ; \
+ } \
+ } \
+}
+
+class BitMaskVec_double {
+
+ MASK_VEC low_, high_ ;
+ SIMD_TYPE combined_ ;
+
+ public:
+ inline MASK_TYPE& getLowEntry(int index) {
+ return low_.masks[index] ;
+ }
+ inline MASK_TYPE& getHighEntry(int index) {
+ return high_.masks[index] ;
+ }
+
+ inline const SIMD_TYPE& getCombinedMask() {
+ VEC_SSE_TO_AVX(low_.vecf, high_.vecf, combined_) ;
+ return combined_ ;
+ }
+
+ inline void shift_left_1bit() {
+ VEC_SHIFT_LEFT_1BIT(low_.vec) ;
+ VEC_SHIFT_LEFT_1BIT(high_.vec) ;
+ }
+
+} ;
+
+#define BITMASK_VEC BitMaskVec_double
diff --git a/public/VectorPairHMM/src/main/c++/define-float.h b/public/VectorPairHMM/src/main/c++/define-float.h
new file mode 100644
index 000000000..318f78280
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/define-float.h
@@ -0,0 +1,206 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#include
+
+#ifdef PRECISION
+#undef PRECISION
+#undef MAIN_TYPE
+#undef MAIN_TYPE_SIZE
+#undef UNION_TYPE
+#undef IF_128
+#undef IF_MAIN_TYPE
+#undef SHIFT_CONST1
+#undef SHIFT_CONST2
+#undef SHIFT_CONST3
+#undef _128_TYPE
+#undef SIMD_TYPE
+#undef AVX_LENGTH
+#undef HAP_TYPE
+#undef MASK_TYPE
+#undef MASK_ALL_ONES
+
+#undef SET_VEC_ZERO(__vec)
+#undef VEC_OR(__v1, __v2)
+#undef VEC_ADD(__v1, __v2)
+#undef VEC_SUB(__v1, __v2)
+#undef VEC_MUL(__v1, __v2)
+#undef VEC_DIV(__v1, __v2)
+#undef VEC_BLEND(__v1, __v2, __mask)
+#undef VEC_BLENDV(__v1, __v2, __maskV)
+#undef VEC_CAST_256_128(__v1)
+#undef VEC_EXTRACT_128(__v1, __im)
+#undef VEC_EXTRACT_UNIT(__v1, __im)
+#undef VEC_SET1_VAL128(__val)
+#undef VEC_MOVE(__v1, __val)
+#undef VEC_CAST_128_256(__v1)
+#undef VEC_INSERT_VAL(__v1, __val, __pos)
+#undef VEC_CVT_128_256(__v1)
+#undef VEC_SET1_VAL(__val)
+#undef VEC_POPCVT_CHAR(__ch)
+#undef VEC_LDPOPCVT_CHAR(__addr)
+#undef VEC_CMP_EQ(__v1, __v2)
+#undef VEC_SET_LSE(__val)
+#undef SHIFT_HAP(__v1, __val)
+#undef MASK_VEC
+#undef VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst)
+#undef VEC_SHIFT_LEFT_1BIT(__vs)
+#undef MASK_ALL_ONES
+#undef COMPARE_VECS(__v1, __v2)
+#undef _256_INT_TYPE
+#undef BITMASK_VEC
+#endif
+
+#define PRECISION s
+
+#define MAIN_TYPE float
+#define MAIN_TYPE_SIZE 32
+#define UNION_TYPE mix_F
+#define IF_128 IF_128f
+#define IF_MAIN_TYPE IF_32
+#define SHIFT_CONST1 12
+#define SHIFT_CONST2 3
+#define SHIFT_CONST3 4
+#define _128_TYPE __m128
+#define SIMD_TYPE __m256
+#define _256_INT_TYPE __m256i
+#define AVX_LENGTH 8
+#define HAP_TYPE UNION_TYPE
+#define MASK_TYPE uint32_t
+#define MASK_ALL_ONES 0xFFFFFFFF
+#define MASK_VEC MaskVec_F
+
+#define SET_VEC_ZERO(__vec) \
+ __vec= _mm256_setzero_ps()
+
+#define VEC_OR(__v1, __v2) \
+ _mm256_or_ps(__v1, __v2)
+
+#define VEC_ADD(__v1, __v2) \
+ _mm256_add_ps(__v1, __v2)
+
+#define VEC_SUB(__v1, __v2) \
+ _mm256_sub_ps(__v1, __v2)
+
+#define VEC_MUL(__v1, __v2) \
+ _mm256_mul_ps(__v1, __v2)
+
+#define VEC_DIV(__v1, __v2) \
+ _mm256_div_ps(__v1, __v2)
+
+#define VEC_BLEND(__v1, __v2, __mask) \
+ _mm256_blend_ps(__v1, __v2, __mask)
+
+#define VEC_BLENDV(__v1, __v2, __maskV) \
+ _mm256_blendv_ps(__v1, __v2, __maskV)
+
+#define VEC_CAST_256_128(__v1) \
+ _mm256_castps256_ps128 (__v1)
+
+#define VEC_EXTRACT_128(__v1, __im) \
+ _mm256_extractf128_ps (__v1, __im)
+
+#define VEC_EXTRACT_UNIT(__v1, __im) \
+ _mm_extract_epi32(__v1, __im)
+
+#define VEC_SET1_VAL128(__val) \
+ _mm_set1_ps(__val)
+
+#define VEC_MOVE(__v1, __val) \
+ _mm_move_ss(__v1, __val)
+
+#define VEC_CAST_128_256(__v1) \
+ _mm256_castps128_ps256(__v1)
+
+#define VEC_INSERT_VAL(__v1, __val, __pos) \
+ _mm256_insertf128_ps(__v1, __val, __pos)
+
+#define VEC_CVT_128_256(__v1) \
+ _mm256_cvtepi32_ps(__v1.i)
+
+#define VEC_SET1_VAL(__val) \
+ _mm256_set1_ps(__val)
+
+#define VEC_POPCVT_CHAR(__ch) \
+ _mm256_cvtepi32_ps(_mm256_set1_epi32(__ch))
+
+#define VEC_LDPOPCVT_CHAR(__addr) \
+ _mm256_cvtepi32_ps(_mm256_loadu_si256((__m256i const *)__addr))
+
+#define VEC_CMP_EQ(__v1, __v2) \
+ _mm256_cmp_ps(__v1, __v2, _CMP_EQ_OQ)
+
+#define VEC_SET_LSE(__val) \
+ _mm256_set_ps(zero, zero, zero, zero, zero, zero, zero, __val);
+
+#define SHIFT_HAP(__v1, __val) \
+ _vector_shift_lastavxs(__v1, __val.f);
+
+#define VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst) \
+ __vdst = _mm256_castps128_ps256(__vsLow) ; \
+__vdst = _mm256_insertf128_ps(__vdst, __vsHigh, 1) ;
+
+#define VEC_SHIFT_LEFT_1BIT(__vs) \
+ __vs = _mm_slli_epi32(__vs, 1)
+
+#define COMPARE_VECS(__v1, __v2, __first, __last) { \
+ float* ptr1 = (float*) (&__v1) ; \
+ float* ptr2 = (float*) (&__v2) ; \
+ for (int ei=__first; ei <= __last; ++ei) { \
+ if (ptr1[ei] != ptr2[ei]) { \
+ std::cout << "Float Mismatch at " << ei << ": " \
+ << ptr1[ei] << " vs. " << ptr2[ei] << std::endl ; \
+ exit(0) ; \
+ } \
+ } \
+}
+
+class BitMaskVec_float {
+
+ MASK_VEC low_, high_ ;
+ SIMD_TYPE combined_ ;
+
+ public:
+
+ inline MASK_TYPE& getLowEntry(int index) {
+ return low_.masks[index] ;
+ }
+ inline MASK_TYPE& getHighEntry(int index) {
+ return high_.masks[index] ;
+ }
+
+ inline const SIMD_TYPE& getCombinedMask() {
+ VEC_SSE_TO_AVX(low_.vecf, high_.vecf, combined_) ;
+ return combined_ ;
+ }
+
+ inline void shift_left_1bit() {
+ VEC_SHIFT_LEFT_1BIT(low_.vec) ;
+ VEC_SHIFT_LEFT_1BIT(high_.vec) ;
+ }
+
+} ;
+
+#define BITMASK_VEC BitMaskVec_float
diff --git a/public/VectorPairHMM/src/main/c++/define-sse-double.h b/public/VectorPairHMM/src/main/c++/define-sse-double.h
new file mode 100644
index 000000000..2d271a854
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/define-sse-double.h
@@ -0,0 +1,173 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#ifdef PRECISION
+#undef PRECISION
+#undef MAIN_TYPE
+#undef MAIN_TYPE_SIZE
+#undef UNION_TYPE
+#undef IF_128
+#undef IF_MAIN_TYPE
+#undef SHIFT_CONST1
+#undef SHIFT_CONST2
+#undef SHIFT_CONST3
+#undef _128_TYPE
+#undef SIMD_TYPE
+#undef AVX_LENGTH
+#undef HAP_TYPE
+#undef MASK_TYPE
+#undef MASK_ALL_ONES
+
+#undef VEC_EXTRACT_UNIT(__v1, __im)
+#undef VEC_INSERT_UNIT(__v1,__ins,__im)
+#undef SET_VEC_ZERO(__vec)
+#undef VEC_OR(__v1, __v2)
+#undef VEC_ADD(__v1, __v2)
+#undef VEC_SUB(__v1, __v2)
+#undef VEC_MUL(__v1, __v2)
+#undef VEC_DIV(__v1, __v2)
+#undef VEC_BLEND(__v1, __v2, __mask)
+#undef VEC_BLENDV(__v1, __v2, __maskV)
+#undef VEC_CAST_256_128(__v1)
+#undef VEC_EXTRACT_128(__v1, __im)
+#undef VEC_EXTRACT_UNIT(__v1, __im)
+#undef VEC_SET1_VAL128(__val)
+#undef VEC_MOVE(__v1, __val)
+#undef VEC_CAST_128_256(__v1)
+#undef VEC_INSERT_VAL(__v1, __val, __pos)
+#undef VEC_CVT_128_256(__v1)
+#undef VEC_SET1_VAL(__val)
+#undef VEC_POPCVT_CHAR(__ch)
+#undef VEC_LDPOPCVT_CHAR(__addr)
+#undef VEC_CMP_EQ(__v1, __v2)
+#undef VEC_SET_LSE(__val)
+#undef SHIFT_HAP(__v1, __val)
+#undef MASK_VEC
+#undef VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst)
+#undef VEC_SHIFT_LEFT_1BIT(__vs)
+#undef MASK_ALL_ONES
+#undef COMPARE_VECS(__v1, __v2)
+#undef _256_INT_TYPE
+#undef BITMASK_VEC
+#endif
+
+#define SSE
+#define PRECISION d
+
+#define MAIN_TYPE double
+#define MAIN_TYPE_SIZE 64
+#define UNION_TYPE mix_D128
+#define IF_128 IF_128d
+#define IF_MAIN_TYPE IF_64
+#define SHIFT_CONST1 1
+#define SHIFT_CONST2 8
+#define SHIFT_CONST3 0
+#define _128_TYPE __m128d
+#define SIMD_TYPE __m128d
+#define _256_INT_TYPE __m128i
+#define AVX_LENGTH 2
+#define HAP_TYPE __m128i
+#define MASK_TYPE uint64_t
+#define MASK_ALL_ONES 0xFFFFFFFFFFFFFFFFL
+#define MASK_VEC MaskVec_D
+
+#define VEC_EXTRACT_UNIT(__v1, __im) \
+ _mm_extract_epi64(__v1, __im)
+
+#define VEC_INSERT_UNIT(__v1,__ins,__im) \
+ _mm_insert_epi64(__v1,__ins,__im)
+
+#define VEC_OR(__v1, __v2) \
+ _mm_or_pd(__v1, __v2)
+
+#define VEC_ADD(__v1, __v2) \
+ _mm_add_pd(__v1, __v2)
+
+#define VEC_SUB(__v1, __v2) \
+ _mm_sub_pd(__v1, __v2)
+
+#define VEC_MUL(__v1, __v2) \
+ _mm_mul_pd(__v1, __v2)
+
+#define VEC_DIV(__v1, __v2) \
+ _mm_div_pd(__v1, __v2)
+
+#define VEC_CMP_EQ(__v1, __v2) \
+ _mm_cmpeq_pd(__v1, __v2)
+
+#define VEC_BLEND(__v1, __v2, __mask) \
+ _mm_blend_pd(__v1, __v2, __mask)
+
+#define VEC_BLENDV(__v1, __v2, __maskV) \
+ _mm_blendv_pd(__v1, __v2, __maskV)
+
+#define SHIFT_HAP(__v1, __val) \
+ __v1 = _mm_insert_epi32(_mm_slli_si128(__v1, 4), __val.i, 0)
+
+#define VEC_CVT_128_256(__v1) \
+ _mm_cvtepi32_pd(__v1)
+
+#define VEC_SET1_VAL(__val) \
+ _mm_set1_pd(__val)
+
+#define VEC_POPCVT_CHAR(__ch) \
+ _mm_cvtepi32_pd(_mm_set1_epi32(__ch))
+
+#define VEC_SET_LSE(__val) \
+ _mm_set_pd(zero, __val);
+
+#define VEC_LDPOPCVT_CHAR(__addr) \
+ _mm_cvtepi32_pd(_mm_loadu_si128((__m128i const *)__addr))
+
+#define VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst) \
+ __vdst = _mm_castsi128_pd(_mm_set_epi64(__vsHigh, __vsLow))
+
+#define VEC_SHIFT_LEFT_1BIT(__vs) \
+ __vs = _mm_slli_epi64(__vs, 1)
+
+
+class BitMaskVec_sse_double {
+
+ MASK_VEC combined_ ;
+ public:
+ inline MASK_TYPE& getLowEntry(int index) {
+ return combined_.masks[index] ;
+ }
+ inline MASK_TYPE& getHighEntry(int index) {
+ return combined_.masks[AVX_LENGTH/2+index] ;
+ }
+
+ inline const SIMD_TYPE& getCombinedMask() {
+ return combined_.vecf ;
+ }
+
+ inline void shift_left_1bit() {
+ VEC_SHIFT_LEFT_1BIT(combined_.vec) ;
+ }
+
+} ;
+
+#define BITMASK_VEC BitMaskVec_sse_double
+
diff --git a/public/VectorPairHMM/src/main/c++/define-sse-float.h b/public/VectorPairHMM/src/main/c++/define-sse-float.h
new file mode 100644
index 000000000..20af947dd
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/define-sse-float.h
@@ -0,0 +1,173 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#ifdef PRECISION
+#undef PRECISION
+#undef MAIN_TYPE
+#undef MAIN_TYPE_SIZE
+#undef UNION_TYPE
+#undef IF_128
+#undef IF_MAIN_TYPE
+#undef SHIFT_CONST1
+#undef SHIFT_CONST2
+#undef SHIFT_CONST3
+#undef _128_TYPE
+#undef SIMD_TYPE
+#undef AVX_LENGTH
+#undef HAP_TYPE
+#undef MASK_TYPE
+#undef MASK_ALL_ONES
+
+#undef VEC_EXTRACT_UNIT(__v1, __im)
+#undef VEC_INSERT_UNIT(__v1,__ins,__im)
+#undef SET_VEC_ZERO(__vec)
+#undef VEC_OR(__v1, __v2)
+#undef VEC_ADD(__v1, __v2)
+#undef VEC_SUB(__v1, __v2)
+#undef VEC_MUL(__v1, __v2)
+#undef VEC_DIV(__v1, __v2)
+#undef VEC_BLEND(__v1, __v2, __mask)
+#undef VEC_BLENDV(__v1, __v2, __maskV)
+#undef VEC_CAST_256_128(__v1)
+#undef VEC_EXTRACT_128(__v1, __im)
+#undef VEC_EXTRACT_UNIT(__v1, __im)
+#undef VEC_SET1_VAL128(__val)
+#undef VEC_MOVE(__v1, __val)
+#undef VEC_CAST_128_256(__v1)
+#undef VEC_INSERT_VAL(__v1, __val, __pos)
+#undef VEC_CVT_128_256(__v1)
+#undef VEC_SET1_VAL(__val)
+#undef VEC_POPCVT_CHAR(__ch)
+#undef VEC_LDPOPCVT_CHAR(__addr)
+#undef VEC_CMP_EQ(__v1, __v2)
+#undef VEC_SET_LSE(__val)
+#undef SHIFT_HAP(__v1, __val)
+#undef MASK_VEC
+#undef VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst)
+#undef VEC_SHIFT_LEFT_1BIT(__vs)
+#undef MASK_ALL_ONES
+#undef COMPARE_VECS(__v1, __v2)
+#undef _256_INT_TYPE
+#undef BITMASK_VEC
+#endif
+
+#define SSE
+#define PRECISION s
+
+#define MAIN_TYPE float
+#define MAIN_TYPE_SIZE 32
+#define UNION_TYPE mix_F128
+#define IF_128 IF_128f
+#define IF_MAIN_TYPE IF_32
+#define SHIFT_CONST1 3
+#define SHIFT_CONST2 4
+#define SHIFT_CONST3 0
+#define _128_TYPE __m128
+#define SIMD_TYPE __m128
+#define _256_INT_TYPE __m128i
+#define AVX_LENGTH 4
+//#define MAVX_COUNT (MROWS+3)/AVX_LENGTH
+#define HAP_TYPE UNION_TYPE
+#define MASK_TYPE uint32_t
+#define MASK_ALL_ONES 0xFFFFFFFF
+#define MASK_VEC MaskVec_F
+
+#define VEC_EXTRACT_UNIT(__v1, __im) \
+ _mm_extract_epi32(__v1, __im)
+
+#define VEC_INSERT_UNIT(__v1,__ins,__im) \
+ _mm_insert_epi32(__v1,__ins,__im)
+
+#define VEC_OR(__v1, __v2) \
+ _mm_or_ps(__v1, __v2)
+
+#define VEC_ADD(__v1, __v2) \
+ _mm_add_ps(__v1, __v2)
+
+#define VEC_SUB(__v1, __v2) \
+ _mm_sub_ps(__v1, __v2)
+
+#define VEC_MUL(__v1, __v2) \
+ _mm_mul_ps(__v1, __v2)
+
+#define VEC_DIV(__v1, __v2) \
+ _mm_div_ps(__v1, __v2)
+
+#define VEC_CMP_EQ(__v1, __v2) \
+ _mm_cmpeq_ps(__v1, __v2)
+
+#define VEC_BLEND(__v1, __v2, __mask) \
+ _mm_blend_ps(__v1, __v2, __mask)
+
+#define VEC_BLENDV(__v1, __v2, __maskV) \
+ _mm_blendv_ps(__v1, __v2, __maskV)
+
+#define SHIFT_HAP(__v1, __val) \
+ _vector_shift_lastsses(__v1, __val.f)
+
+#define VEC_CVT_128_256(__v1) \
+ _mm_cvtepi32_ps(__v1.i)
+
+#define VEC_SET1_VAL(__val) \
+ _mm_set1_ps(__val)
+
+#define VEC_POPCVT_CHAR(__ch) \
+ _mm_cvtepi32_ps(_mm_set1_epi32(__ch))
+
+#define VEC_SET_LSE(__val) \
+ _mm_set_ps(zero, zero, zero, __val);
+
+#define VEC_LDPOPCVT_CHAR(__addr) \
+ _mm_cvtepi32_ps(_mm_loadu_si128((__m128i const *)__addr))
+
+#define VEC_SSE_TO_AVX(__vsLow, __vsHigh, __vdst) \
+ __vdst = _mm_cvtpi32x2_ps(__vsLow, __vsHigh)
+
+#define VEC_SHIFT_LEFT_1BIT(__vs) \
+ __vs = _mm_slli_epi32(__vs, 1)
+
+class BitMaskVec_sse_float {
+
+ MASK_VEC combined_ ;
+
+ public:
+ inline MASK_TYPE& getLowEntry(int index) {
+ return combined_.masks[index] ;
+ }
+ inline MASK_TYPE& getHighEntry(int index) {
+ return combined_.masks[AVX_LENGTH/2+index] ;
+ }
+
+ inline const SIMD_TYPE& getCombinedMask() {
+ return combined_.vecf ;
+ }
+
+ inline void shift_left_1bit() {
+ VEC_SHIFT_LEFT_1BIT(combined_.vec) ;
+ }
+
+} ;
+
+#define BITMASK_VEC BitMaskVec_sse_float
diff --git a/public/VectorPairHMM/src/main/c++/headers.h b/public/VectorPairHMM/src/main/c++/headers.h
new file mode 100644
index 000000000..4a0d89b57
--- /dev/null
+++ b/public/VectorPairHMM/src/main/c++/headers.h
@@ -0,0 +1,71 @@
+/*Copyright (c) 2012 The Broad Institute
+
+*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
+*NONINFRINGEMENT. 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.
+*/
+
+
+#ifndef COMMON_HEADERS_H
+#define COMMON_HEADERS_H
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include