beagleutil/ChromInterval.java

@@ -205,7 +205,7 @@ public final class ChromInterval implements IntInterval,
      * @return the last genome coordinate in this chromosome interval.
      */
     @Override
-    public int end() {
+    public int inclEnd() {
         return end;
     }
 
@@ -332,7 +332,7 @@ public final class ChromInterval implements IntInterval,
             return false;
         }
         else {
-            return (a.start() <= b.end()) && (b.start() <= a.end());
+            return (a.start() <= b.inclEnd()) && (b.start() <= a.inclEnd());
         }
     }
 
@@ -350,7 +350,7 @@ public final class ChromInterval implements IntInterval,
             throw new IllegalArgumentException(s);
         }
         int start = Math.min(a.start(), b.start());
-        int end = Math.max(a.end(), b.end());
+        int end = Math.max(a.inclEnd(), b.inclEnd());
         return new ChromInterval(a.chrom(), start, end);
     }
 }

beagleutil/IntInterval.java

@@ -40,13 +40,13 @@ public interface IntInterval {
      * Returns the end of the interval (inclusive).
      * @return the end of the interval (inclusive).
      */
-    public int end();
+    public int inclEnd();
 
     /**
      * Returns a {@code Comparator<IntInterval>} which orders
      * {@code IntInterval} objects in order of increasing {@code this.start()}
      * value and orders {@code IntInterval} objects with the same
-     * {@code this.start()} value in order of increasing {@code this.end()}
+     * {@code this.start()} value in order of increasing {@code this.inclEnd()}
      * value.
      * @return a {@code Comparator<IntInterval>} object
      */
@@ -55,8 +55,8 @@ public interface IntInterval {
             if (t1.start() != t2.start()) {
                 return (t1.start() < t2.start()) ? -1 : 1;
             }
-            else if (t1.end() != t2.end()) {
-                return (t1.end() < t2.end()) ? -1 : 1;
+            else if (t1.inclEnd() != t2.inclEnd()) {
+                return (t1.inclEnd() < t2.inclEnd()) ? -1 : 1;
             }
             return 0;
         } ;
@@ -66,7 +66,7 @@ public interface IntInterval {
      * Returns a {@code Comparator<IntInterval>} which orders
      * {@code IntInterval} objects in order of increasing {@code this.start()}
      * value and orders {@code IntInterval} objects with the same
-     * {@code this.start()} value in order of decreasing {@code this.end()}
+     * {@code this.start()} value in order of decreasing {@code this.inclEnd()}
      * value.
      * @return a {@code Comparator<IntInterval>} object
      */
@@ -75,8 +75,8 @@ public interface IntInterval {
             if (t1.start() != t2.start()) {
                 return (t1.start() < t2.start()) ? -1 : 1;
             }
-            else if (t1.end() != t2.end()) {
-                return (t1.end() > t2.end()) ? -1 : 1;
+            else if (t1.inclEnd() != t2.inclEnd()) {
+                return (t1.inclEnd() > t2.inclEnd()) ? -1 : 1;
             }
             return 0;
         } ;

beagleutil/PbwtUpdater.java

+/*
+ * Copyright (C) 2014-2016 Brian L. Browning
+ *
+ * This file is part of Beagle
+ *
+ * Beagle is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Beagle is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+package beagleutil;
+
+import ints.IntArray;
+import ints.IntList;
+import java.util.Arrays;
+import java.util.stream.IntStream;
+
+/**
+ * <p>Class {@code PbwtUpdater} updates prefix and divergence arrays using
+ * the positional Burrows-Wheeler transform (PBWT).</p>
+ *
+ * <p>Instances of {@code PbwtUpdater} are not thread-safe.</p>
+ *
+ * <p>Reference: Durbin, Richard (2014) Efficient haplotype matching and storage
+ *    using the positional Burrows-Wheeler transform (PBWT).
+ *    Bioinformatics 30(9):166-1272. doi: 10.1093/bioinformatics/btu014</p>
+ *
+ * @author Brian L. Browning {@code <browning@uw.edu>}
+ */
+public class PbwtUpdater {
+
+    private final int nHaps;
+
+    private IntList[] a; // next prefix array data
+    private IntList[] d; // next div array data
+    private int[] p;
+
+    /**
+     * Constructs a new {@code PbwtUpdater} instance for the specified data.
+     * @param nHaps the number of haplotypes at each position
+     * @throws NegativeArraySizeException if {@code nHaps < 0}
+     */
+    public PbwtUpdater(int nHaps) {
+        int initSize = 4;
+        this.nHaps = nHaps;
+        this.p = new int[initSize];
+        this.a = IntStream.range(0, initSize)
+                .mapToObj(i -> new IntList())
+                .toArray(IntList[]::new);
+        this.d = IntStream.range(0, initSize)
+                .mapToObj(i -> new IntList())
+                .toArray(IntList[]::new);
+    }
+
+    /**
+     * Returns the number of haplotypes.
+     * @return the number of haplotypes
+     */
+    public int nHaps() {
+        return nHaps;
+    }
+
+    /**
+     * Update the specified prefix and divergence arrays using the forward
+     * Positional Burrows-Wheeler Transform. The contract for this method is
+     * undefined if the specified {@code prefix} array is not a permutation of
+     * {@code 0, 1, 2, ..., (nHaps - 1)}, or if the elements of the
+     * specified {@code div} arrays are not all less than or equal to
+     * the specified {@code marker}.
+     *
+     * @param rec the haplotype alleles
+     * @param nAlleles the number of alleles
+     * @param marker the marker index for the specified haplotype alleles
+     * @param prefix the prefix array
+     * @param div the divergence array
+     *
+     * @throws IllegalArgumentException if {@code nAlleles < 1}
+     * @throws IndexOutOfBoundsException if
+     * {@code (rec.get[j] < 0 || rec.get(j) >= nAlleles)}
+     * for any {@code j} satisfying {@code (0 <= j && j < this.nHaps())}
+     * @throws IndexOutOfBoundsException if {@code prefix.length >= this.nHaps()}
+     * @throws IndexOutOfBoundsException if {@code div.length >= this.nHaps()}
+     * @throws NullPointerException if
+     * {@code rec == null || prefix == null || div == null}
+     */
+    public void fwdUpdate(IntArray rec, int nAlleles, int marker, int[] prefix,
+            int[] div) {
+        initializeArrays(nAlleles, marker+1);
+        for (int i=0; i<nHaps; ++i) {
+            int allele = rec.get(prefix[i]);
+            if (allele>=nAlleles) {
+                throw new IndexOutOfBoundsException(String.valueOf(nAlleles));
+            }
+            for (int j=0; j<nAlleles; ++j) {
+                if (div[i]>p[j]) {
+                    p[j] = div[i];
+                }
+            }
+            a[allele].add(prefix[i]);
+            d[allele].add(p[allele]);
+            p[allele] = Integer.MIN_VALUE;
+        }
+        updatePrefixAndDiv(nAlleles, prefix, div);
+    }
+
+    /**
+     * Update the specified prefix and divergence arrays using the backward
+     * Positional Burrows-Wheeler Transform. The contract for this method is
+     * undefined if the specified {@code prefix} array is not a permutation of
+     * {@code 0, 1, 2, ..., (nHaps - 1)}, or if the elements of the
+     * specified {@code div} arrays are not all greater than or equal to
+     * the specified {@code marker}.
+     *
+     * @param rec the haplotype alleles
+     * @param nAlleles the number of alleles
+     * @param marker the marker index for the specified haplotype alleles
+     * @param prefix the prefix array
+     * @param div the divergence array
+     *
+     * @throws IllegalArgumentException if {@code nAlleles < 1}
+     * @throws IndexOutOfBoundsException if
+     * {@code (rec.get[j] < 0 || rec.get(j) >= nAlleles)}
+     * for any {@code j} satisfying {@code (0 <= j && j < this.nHaps())}
+     * @throws IndexOutOfBoundsException if {@code prefix.length >= this.nHaps()}
+     * @throws IndexOutOfBoundsException if {@code div.length >= this.nHaps()}
+     * @throws NullPointerException if
+     * {@code rec == null || prefix == null || div == null}
+     */
+    public void bwdUpdate(IntArray rec, int nAlleles, int marker, int[] prefix,
+            int[] div) {
+        initializeArrays(nAlleles, marker-1);
+        for (int i=0; i<nHaps; ++i) {
+            int allele = rec.get(prefix[i]);
+            if (allele>=nAlleles) {
+                throw new IndexOutOfBoundsException(String.valueOf(nAlleles));
+            }
+            for (int j=0; j<nAlleles; ++j) {
+                if (div[i]<p[j]) {
+                    p[j] = div[i];
+                }
+            }
+            a[allele].add(prefix[i]);
+            d[allele].add(p[allele]);
+            p[allele] = Integer.MAX_VALUE;
+        }
+        updatePrefixAndDiv(nAlleles, prefix, div);
+    }
+
+    private void updatePrefixAndDiv(int nAlleles, int[] prefix, int[] div) {
+        int start = 0;
+        for (int al=0; al<nAlleles; ++al) {
+            int size = a[al].size();
+            System.arraycopy(a[al].toArray(), 0, prefix, start, size);
+            System.arraycopy(d[al].toArray(), 0, div, start, size);
+            start += size;
+            a[al].clear();
+            d[al].clear();
+        }
+        assert start == nHaps;
+    }
+
+    private void initializeArrays(int nAlleles, int initPValue) {
+        if (nAlleles<1) {
+            throw new IllegalArgumentException(String.valueOf(nAlleles));
+        }
+        ensureArrayCapacity(nAlleles);
+        Arrays.fill(p, 0, nAlleles, initPValue);
+    }
+
+    private void ensureArrayCapacity(int nAlleles) {
+        if (nAlleles>a.length) {
+            int oldLength = a.length;
+            p = Arrays.copyOf(p, nAlleles);
+            a = Arrays.copyOf(a, nAlleles);
+            d = Arrays.copyOf(d, nAlleles);
+            for (int j = oldLength; j<a.length; ++j) {
+                a[j] = new IntList();
+                d[j] = new IntList();
+            }
+        }
+    }
+}

haplotype/HapPair.java → blbutil/DoubleArray.java

@@ -16,85 +16,77 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-package haplotype;
+package blbutil;
 
-import beagleutil.Samples;
-import java.util.Comparator;
-import vcf.Marker;
-import vcf.Markers;
+import java.util.Arrays;
 
 /**
- * <p>Interface {@code HapPair} represents a pair of haplotypes for a sample.
- * The pair of haplotypes are guaranteed to have non-missing alleles at each
- * marker.
- * </p>
- * All instances of {@code HapPair} are required to be immutable.
+ * Class {@code DoubleArray} represents an immutable list of double floating
+ * point values.
  *
  * @author Brian L. Browning {@code <browning@uw.edu>}
  */
-public interface HapPair {
+public class DoubleArray {
 
-    /**
-     * Returns the first allele for the specified marker.
-     * @param marker a marker index
-     * @return the first allele for the specified marker
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
-     */
-    int allele1(int marker);
+    private final double[] values;
 
     /**
-     * Returns the second allele for the specified marker.
-     * @param marker a marker index
-     * @return the second allele for the specified marker
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
+     * Constructs an {@code DoubleArray} object with the specified
+     * values.
+     * @param values the list of floating point values
+     * @throws NullPointerException if {@code values == null}
      */
-    int allele2(int marker);
+    public DoubleArray(double[] values) {
+        this.values = values.clone();
+    }
 
     /**
-     * Returns the markers.
-     * @return the markers
+     * Returns the double at the specified position in this list.
+     * @param index the index of the returned double
+     * @return the double at the specified position in this list
+     * @throws IndexOutOfBoundsException if
+     * {@code index < 0 || index >= this.size()}
      */
-    Markers markers();
+    public double get(int index) {
+        return values[index];
+    }
 
     /**
-     * Returns the specified marker.
-     * @param marker a marker index
-     * @return the specified marker
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
+     * Returns the number of elements in this list.
+     * @return the number of elements in this list
      */
-     Marker marker(int marker);
+    public int size() {
+        return values.length;
+    }
 
     /**
-     * Returns the number of markers.
-     * @return the number of markers
+     * Returns {@code true} if this list has no elements, and returns
+     * {@code false} otherwise.
+     * @return {@code true} if this list has no elements, and returns
+     * {@code false} otherwise
      */
-    int nMarkers();
+    public boolean isEmpty() {
+        return values.length==0;
+    }
 
     /**
-     * Returns the sample identifier index.
-     * @return the sample identifier index
+     * Returns an integer array containing the sequence of elements in this
+     * list.
+     * @return an integer array containing the sequence of elements in this
+     * list
      */
-    int idIndex();
+    public double[] toArray() {
+        return values.clone();
+    }
 
     /**
-     * Returns a {@code Comparator<HapPairInterface>}
-     * whose {@code compare(hp1, hp2)} method returns -1, 0, or 1
-     * depending on whether {@code samples.index(hp1.idIndex())} is
-     * less than, equal, or greater than
-     * {@code samples.index(hp2.idIndex())}.
-     * @param samples the list of samples used to compare {@code HapsPair}
-     * objects
-     * @return a {@code Comparator<HapPairInterface>}
-     * whose {@code compare(hp1, hp2)} method compares two
-     * haplotype pairs for order
-     * @throws NullPointerException if {@code samples == null}
+     * Returns a string representation of this list that is
+     * obtained by calling {@code java.util.Arrays.toString(this.toArray())}.
+     *
+     * @return a string representation of this list
      */
-    static Comparator<HapPair> comparator(final Samples samples) {
-        return (hp1, hp2) -> Integer.compare(
-                samples.index(hp1.idIndex()),
-                samples.index(hp2.idIndex()));
+    @Override
+    public String toString() {
+        return Arrays.toString(values);
     }
 }

blbutil/FloatArray.java

@@ -40,6 +40,19 @@ public class FloatArray {
         this.values = values.clone();
     }
 
+    /**
+     * Constructs an {@code FloatArray} object with the specified
+     * values.
+     * @param values the list of floating point values
+     * @throws NullPointerException if {@code values == null}
+     */
+    public FloatArray(double[] values) {
+        this.values = new float[values.length];
+        for (int j=0; j<values.length; ++j) {
+            this.values[j] = (float) values[j];
+        }
+    }
+
     /**
      * Returns the float at the specified position in this list.
      * @param index the index of the returned float

blbutil/Utilities.java

@@ -249,10 +249,8 @@ public class Utilities {
      * @throws NullPointerException if {@code e == null}
      */
     public static void exit(Throwable t) {
-        System.out.print("ERROR: ");
-        System.out.println(t);
+        System.out.println("ERROR: " + t);
         t.printStackTrace(System.out);
-        System.out.println("terminating program.");
         System.exit(1);
     }
 

bref/Bref3Reader.java

@@ -35,7 +35,7 @@ import java.util.List;
 import vcf.BasicMarker;
 import vcf.Marker;
 import vcf.RefGTRec;
-import vcf.SeqCodedRefGT;
+import vcf.SeqCodedRefGTRec;
 
 /**
  * <p>Class {@code Bref3Reader} contains methods for reading a bref 3
@@ -180,7 +180,7 @@ public final class Bref3Reader {
 //            throw new IllegalArgumentException("inconsistent data");
 //        }
 
-        return new SeqCodedRefGT(marker, samples, hapToSeq, seqToAllele);
+        return new SeqCodedRefGTRec(marker, samples, hapToSeq, seqToAllele);
     }
 
     private static void readAndCheckMagicNumber(DataInput di) throws IOException {

bref/SeqCoder3.java

@@ -27,7 +27,7 @@ import java.util.Collections;
 import java.util.List;
 import vcf.Marker;
 import vcf.RefGTRec;
-import vcf.SeqCodedRefGT;
+import vcf.SeqCodedRefGTRec;
 
 /**
  * <p>Class {@code SeqCoder3} compresses a sequence of allele-coded
@@ -286,7 +286,7 @@ public class SeqCoder3 {
             RefGTRec rec = recs.get(j);
             Marker m = rec.marker();
             IntArray seq2allele = seq2Allele(rec, seq2Hap);
-            list.add(new SeqCodedRefGT(m, samples, hap2seq, seq2allele));
+            list.add(new SeqCodedRefGTRec(m, samples, hap2seq, seq2allele));
         }
         initialize();
         return list;

debian/beagle-doc.doc-base

@@ -3,7 +3,7 @@ Title: Beagle
 Section: Science/Biology
 
 Format: PDF
-Files: /usr/share/doc/beagle/beagle_5.0.pdf.gz
+Files: /usr/share/doc/beagle/beagle_5.1.pdf.gz
 
 Format: HTML
 Index: /usr/share/doc/beagle/api/index.html

debian/beagle-doc.install

 debian/tests/run-sample-analysis /usr/share/doc/beagle/examples/
 debian/upstream.docs/run.beagle.example /usr/share/doc/beagle/examples/
 debian/upstream.docs/test.vcf /usr/share/doc/beagle/examples/
-debian/upstream.docs/beagle_5.0.pdf /usr/share/doc/beagle/
+debian/upstream.docs/beagle_5.1.pdf /usr/share/doc/beagle/
 debian/_jh_build.javadoc/api /usr/share/doc/beagle

debian/changelog

+beagle (5.1-190824+dfsg-1) UNRELEASED; urgency=medium
+
+  * New upstream release.
+
+ -- Dylan Aïssi <daissi@debian.org>  Mon, 09 Sep 2019 21:58:35 +0200
+
 beagle (5.0-190712+dfsg-1) unstable; urgency=medium
 
   * New upstream release.

debian/patches/01_Update_htsjtk.patch

-Author: Dylan Aïssi <bob.dybian@gmail.com>
+Author: Dylan Aïssi <daissi@debian.org>
 Description: Update references to htsjdk.
-Last-Update: 2016-04-09
+Last-Update: 2019-09-09
 Forwarded: No
 
 --- a/blbutil/InputIt.java

debian/rules

@@ -8,7 +8,7 @@ export CLASSPATH=/usr/share/java/htsjdk.jar
 
 override_jh_build:
 	mkdir SRC/
-	cp -t SRC/ -r beagleutil blbutil bref haplotype imp ints main math phase sample vcf
+	cp -t SRC/ -r beagleutil blbutil bref imp ints main math phase vcf
 	jh_build --javacopts="-source 1.8" --javadoc-opts="-source 1.8" --main="main.Main" beagle.jar SRC/
 	jh_build --javacopts="-source 1.8" --javadoc-opts="-source 1.8" --main="bref.Bref3" bref3.jar SRC/
 

debian/source/include-binaries

 debian/upstream.docs/beagle_4.1_09Feb16.docx
-debian/upstream.docs/beagle_4.1.pdf
-debian/upstream.docs/beagle_5.0.pdf
+debian/upstream.docs/beagle_5.1.pdf

debian/upstream.docs/README

 # Author: Dylan Assi <daissi@debian.org>
-# Last-Update: 2019-07-13
+# Last-Update: 2019-09-09
 # These files were downloaded from https://faculty.washington.edu/browning/beagle/
 # They should be updated when upstream update them.
 
 	wget https://faculty.washington.edu/browning/beagle/beagle5_release_notes -O release_notes
-	wget https://faculty.washington.edu/browning/beagle/beagle_5.0_03Jul19.pdf -O beagle_5.0.pdf
-	wget https://faculty.washington.edu/browning/beagle/run.beagle.12Jul19.0df.example -O run.beagle.example
-	wget http://faculty.washington.edu/browning/beagle/test.12Jul19.0df.vcf.gz -O test.vcf.gz && \
+	wget https://faculty.washington.edu/browning/beagle/beagle_5.1_12Aug19.pdf -O beagle_5.1.pdf
+	wget https://faculty.washington.edu/browning/beagle/run.beagle.24Aug19.3e8.example -O run.beagle.example
+	wget https://faculty.washington.edu/browning/beagle/test.24Aug19.3e8.vcf.gz -O test.vcf.gz && \
 		gunzip test.vcf.gz

debian/upstream.docs/release_notes

@@ -64,3 +64,12 @@ Beagle 5.0 (12Jul19.0df) release notes
 =============================================
 * Improved an error message
 * Improved handling of inconsistent marker order in ref and gt files
+
+Beagle 5.1 (13Aug19.a31) release notes
+============================================
+* Initial release of Beagle 5.1
+* Beagle 5.1 has some improvements that increase speed and accuracy
+
+Beagle 5.1 (24Aug19.3e8) release notes
+============================================
+* Fixed a bug that can arise when phasing inbred samples

debian/upstream.docs/run.beagle.example

 #!/bin/bash
-if [ ! -f beagle.12Jul19.0df.jar ]; then
+if [ ! -f beagle.24Aug19.3e8.jar ]; then
   echo
-  echo "Downloading beagle.12Jul19.0df.jar"
-  wget http://faculty.washington.edu/browning/beagle/beagle.12Jul19.0df.jar
+  echo "Downloading beagle.24Aug19.3e8.jar"
+  wget http://faculty.washington.edu/browning/beagle/beagle.24Aug19.3e8.jar
 fi
 
-if [ ! -f bref3.12Jul19.0df.jar ]; then
+if [ ! -f bref3.24Aug19.3e8.jar ]; then
   echo
-  echo "Downloading bref3.12Jul19.0df.jar"
-  wget http://faculty.washington.edu/browning/beagle/bref3.12Jul19.0df.jar
+  echo "Downloading bref3.24Aug19.3e8.jar"
+  wget http://faculty.washington.edu/browning/beagle/bref3.24Aug19.3e8.jar
 fi
 
 echo
 
-if [ ! -f test.12Jul19.0df.vcf.gz ]; then
+if [ ! -f test.24Aug19.3e8.vcf.gz ]; then
     echo
-    echo "*** Downloading some 1000 Genomes Project data to file: test.12Jul19.0df.vcf.gz ***"
-    wget http://faculty.washington.edu/browning/beagle/test.12Jul19.0df.vcf.gz
+    echo "*** Downloading some 1000 Genomes Project data to file: test.24Aug19.3e8.vcf.gz ***"
+    wget http://faculty.washington.edu/browning/beagle/test.24Aug19.3e8.vcf.gz
 fi
 
 echo
-echo "*** Creating test files: ref.12Jul19.0df.vcf.gz target.12Jul19.0df.vcf.gz ***"
+echo "*** Creating test files: ref.24Aug19.3e8.vcf.gz target.24Aug19.3e8.vcf.gz ***"
 echo
-zcat test.12Jul19.0df.vcf.gz | cut -f1-190 | tr '/' '|' | gzip > ref.12Jul19.0df.vcf.gz
-zcat test.12Jul19.0df.vcf.gz | cut -f1-9,191-200 | gzip > target.12Jul19.0df.vcf.gz
+zcat test.24Aug19.3e8.vcf.gz | cut -f1-190 | tr '/' '|' | gzip > ref.24Aug19.3e8.vcf.gz
+zcat test.24Aug19.3e8.vcf.gz | cut -f1-9,191-200 | gzip > target.24Aug19.3e8.vcf.gz
 
 echo
 echo "*** Running test analysis with \"gt=\" argument ***"
 echo
-java -jar beagle.12Jul19.0df.jar gt=test.12Jul19.0df.vcf.gz out=out.gt
+java -jar beagle.24Aug19.3e8.jar gt=test.24Aug19.3e8.vcf.gz out=out.gt
 
 echo
 echo "*** Running test analysis with \"ref=\" and \"gt=\" arguments ***"
 echo
-java -jar beagle.12Jul19.0df.jar ref=ref.12Jul19.0df.vcf.gz gt=target.12Jul19.0df.vcf.gz out=out.ref
+java -jar beagle.24Aug19.3e8.jar ref=ref.24Aug19.3e8.vcf.gz gt=target.24Aug19.3e8.vcf.gz out=out.ref
 
 echo
 echo "*** Making \"bref3\" file ***"
 echo
-java -jar bref3.12Jul19.0df.jar ref.12Jul19.0df.vcf.gz > ref.12Jul19.0df.bref3
+java -jar bref3.24Aug19.3e8.jar ref.24Aug19.3e8.vcf.gz > ref.24Aug19.3e8.bref3
 
 echo
 echo "*** Running test analysis with \"bref3\" file ***"
 echo
-java -jar beagle.12Jul19.0df.jar ref=ref.12Jul19.0df.bref3 gt=target.12Jul19.0df.vcf.gz out=out.bref3
+java -jar beagle.24Aug19.3e8.jar ref=ref.24Aug19.3e8.bref3 gt=target.24Aug19.3e8.vcf.gz out=out.bref3
 

haplotype/BitHapPair.java

-/*
- * Copyright (C) 2014-2016 Brian L. Browning
- *
- * This file is part of Beagle
- *
- * Beagle is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Beagle is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-package haplotype;
-
-import vcf.Marker;
-import vcf.Markers;
-
-/**
- * <p>Class {@code BitHapPair} represents a pair of haplotypes for a sample.
- * The class stores alleles as a bit array.
- * </p>
- * Instances of class {@code BitHapPair} are immutable.
- *
- * @author Brian L. Browning {@code <browning@uw.edu>}
- */
-public final class BitHapPair implements HapPair {
-
-    private static final int LOG2_BITS_PER_WORD = 6;
-
-    private final Markers markers;
-    private final int idIndex;
-    private final long[] bits1;
-    private final long[] bits2;
-
-    /**
-     * Constructs a new {@code BitHapPair} instance.
-     * @param markers the sequence of markers
-     * @param idIndex the sample identifier index
-     * @param alleles1 the sequence of allele indices for the first haplotype
-     * @param alleles2 the sequence of alleles indices for the second haplotype
-     *
-     * @throws IllegalArgumentException if
-     * {@code alleles1.length != markers.nMarkers()
-     * || alleles2.length != markers.nMarkers()}
-     * @throws IllegalArgumentException if {@code alleles1[k] < 0 ||
-     * allele1[k] >= markers.marker(k).nAlleles()} for some {@code k} satisfying
-     * {@code 0 <= k && k < markers.nMarkers()}
-     * @throws IllegalArgumentException if {@code alleles2[k] < 0 ||
-     * allele2[k] >= markers.marker(k).nAlleles()} for some {@code k} satisfying
-     * {@code 0 <= k && k < markers.nMarkers()}
-     * @throws IndexOutOfBoundsException if {@code idIndex < 0}
-     * @throws NullPointerException if
-     * {@code marker == null || alleles1 == null || allele2 == null}
-     */
-    public BitHapPair(Markers markers, int idIndex, int[] alleles1,
-            int[] alleles2) {
-        if (alleles1.length != markers.nMarkers()
-                || alleles2.length != markers.nMarkers()) {
-            throw new IllegalArgumentException("inconsistent markers");
-        }
-        if (idIndex < 0) {
-            throw new IndexOutOfBoundsException(String.valueOf(idIndex));
-        }
-        this.markers = markers;
-        this.idIndex = idIndex;
-        this.bits1 = toBitArray(markers, alleles1);
-        this.bits2 = toBitArray(markers, alleles2);
-    }
-
-    private static long[] toBitArray(Markers markers, int[] alleles) {
-        int nWords = (markers.sumHaplotypeBits() + (Long.SIZE-1)) >> LOG2_BITS_PER_WORD;
-        long[] bits = new long[nWords];
-        int bitIndex = 0;
-        for (int k=0; k<alleles.length; ++k) {
-            int allele = alleles[k];
-            if (allele < 0 || allele >= markers.marker(k).nAlleles()) {
-                String s = "allele \"" + allele + "\" out of bounds for marker: "
-                        + markers.marker(k);
-                throw new IllegalArgumentException(s);
-            }
-            int mask = 1;
-            int nBits = markers.sumHaplotypeBits(k+1) - markers.sumHaplotypeBits(k);
-            for (int l=0; l<nBits; ++l) {
-                if ((allele & mask)==mask) {
-                    int wordIndex =  bitIndex >> LOG2_BITS_PER_WORD;
-                    bits[wordIndex] |= (1L << bitIndex);
-                }
-                bitIndex++;
-                mask <<= 1;
-            }
-        }
-        return bits;
-    }
-
-    @Override
-    public int allele1(int marker) {
-        return allele(bits1, marker);
-    }
-
-    @Override
-    public int allele2(int marker) {
-        return allele(bits2, marker);
-    }
-
-    private int allele(long[] alleleBits, int marker) {
-        int start = markers.sumHaplotypeBits(marker);
-        int end = markers.sumHaplotypeBits(marker+1);
-        if (end==(start+1)) {
-            int wordIndex =  start >> LOG2_BITS_PER_WORD;
-            return (int) (alleleBits[wordIndex] >> start) & 1;
-        }
-        int allele = 0;
-        int mask = 1;
-        for (int j=start; j<end; ++j) {
-            int wordIndex =  j >> LOG2_BITS_PER_WORD;
-            if ((alleleBits[wordIndex] & (1L << j)) != 0) {
-                allele |= mask;
-            }
-            mask <<= 1;
-        }
-        return allele;
-    }
-
-    @Override
-    public Markers markers() {
-        return markers;
-    }
-
-    @Override
-    public Marker marker(int marker) {
-        return markers.marker(marker);
-    }
-
-    @Override
-    public int nMarkers() {
-        return markers.nMarkers();
-    }
-
-    @Override
-    public int idIndex() {
-        return idIndex;
-    }
-
-    /**
-     * Returns a string representation of {@code this}.  The
-     * exact details of the representation are unspecified and subject
-     * to change.
-     * @return a string representation of {@code this}
-     */
-    @Override
-    public String toString() {
-        StringBuilder sb = new StringBuilder();
-        sb.append("BitHapPair: idIndex=");
-        sb.append(idIndex);
-        return sb.toString();
-    }
-}