beagleutil/Phase.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;
-
-/**
- * <p>Class {@code Phase} represents the equivalence of two phased genotypes
- * for a marker or for a set of markers.  Genotype equivalence is defined
- * in terms of allele equivalence.  Two alleles are equivalent if
- * either allele is missing or if both alleles are non-missing and equal.
- * </p>
- * <p>
- * For the case of a single marker with phased (i.e. ordered) genotypes
- * ({@code a1}, {@code a2})
- * and ({@code b1}, {@code b2}), then
- * <br>
- * 1) the genotypes have IDENTICAL phase if a) alleles {@code a1}
- * and {@code b1} are equivalent, b) alleles {@code a2} and
- * {@code b2} are equivalent, and c) either alleles {@code a1}
- * and {@code b2} are not equivalent or alleles
- * {@code a2} and {@code b1} are not equivalent.
- * <br>
- * 2) the genotypes have OPPOSITE phase if a) alleles {@code a1}
- * and {@code b2} are equivalent, b) alleles {@code a2} and
- * {@code b1} are equivalent, and c) either alleles {@code a1}
- * and {@code b1} are not equivalent or alleles {@code a2} and
- * {@code b2} are not equivalent.
- * <br>
- * 3) the genotypes have UNKOWN phase if a) alleles {@code a1}
- * and {@code b1} are equivalent, b) alleles {@code a2} and
- * {@code b2} are equivalent, c) alleles {@code a1} and
- * {@code b2} are equivalent, and d) alleles {@code a2} and
- * {@code b1} are equivalent.
- * <br>
- * 4) the genotypes have INCONSISTENT phase if a) either alleles
- * {@code a1} and {@code b1} are not equivalent or alleles
- * {@code a2} and {@code b2} are not equivalent, and
- * b) either alleles {@code a1} and {@code b2} are not equivalent
- * or alleles {@code a2} and {@code b1} are not equivalent.
- * </p>
- * For the case of two sets of phased genotypes for the same markers,
- * the two sets have
- * <br>
- * 1) IDENTICAL phase if the phase is IDENTICAL for at least
- * one marker and is either IDENTICAL or UNKNOWN for all markers.
- * <br>
- * 2) OPPOSITE phase if the if the phase is OPPOSITE for at
- * lease one marker and is either OPPOSITE or UNKNOWN for all markers.
- * <br>
- * 3) UNKNOWN phase if the phase is UNKNOWN for all markers.
- * <br>
- * 4) INCONSISTENT phase if a) the phase is INCONSISTENT for at least one
- * marker or if b) the relative phase is IDENTICAL for at least one marker and
- * OPPOSITE for at least one marker.
- *
- * @author Brian L. Browning
- */
-public enum Phase {
-        IDENTICAL,
-        OPPOSITE,
-        UNKNOWN,
-        INCONSISTENT
-}

blbutil/MultiThreadUtils.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 blbutil;
+
+import java.util.concurrent.BlockingQueue;
+import java.util.concurrent.CountDownLatch;
+import java.util.concurrent.ExecutorService;
+import java.util.concurrent.TimeUnit;
+
+/**
+ * Class {@code Utilities} contains miscellaneous static utility methods
+ * for multi-threaded programming.
+ *
+ * @author Brian L. Browning {@code <browning@uw.edu>}
+ */
+public class MultiThreadUtils {
+
+    private MultiThreadUtils() {
+        // private constructor to prevent instantiation
+    }
+
+    /**
+     * Inserts the specified element at the tail of the specified blocking
+     * queue, waiting for space to become available if the queue is full.
+     * The Java Virtual Machine is terminated if an {@code InterruptedException}
+     * is thrown while waiting for space to be come available in the queue.
+     * @param <E> the element type
+     * @param q a blocking queue
+     * @param e the element to add
+     * @throws NullPointerException if {@code q == null || e == null}
+     */
+    public static <E> void putInBlockingQ(BlockingQueue<E> q, E e) {
+        try {
+            q.put(e);
+        } catch (InterruptedException ex) {
+            Utilities.exit("ERROR: " , ex);
+        }
+    }
+
+    /**
+     * Removes and returns the element at the head of the specified blocking
+     * queue, waiting if necessary for an element to become available.
+     * The Java Virtual Machine is terminated if an {@code InterruptedException}
+     * is thrown while waiting for space to be come available in the queue.
+     * @param <E> the element type
+     * @param q a blocking queue
+     * @return the element at the head of the queue
+     */
+    public static <E> E takeFromBlockingQ(BlockingQueue<E> q) {
+        try {
+            return q.take();
+        } catch (InterruptedException ex) {
+            Utilities.exit("ERROR: " , ex);
+        }
+        assert false;
+        return null;
+    }
+
+    /**
+     * Blocks the current thread until the specified {@code CountDownLatch}
+     * has counted down to 0. The Java Virtual Machine is terminated if an
+     * {@code InterruptedException} is thrown while waiting for for the
+     * {@code CountDownLatch} to count down to 0.
+     * @param latch the count down latch
+     * @throws NullPointerException if {@code latch == null}
+     */
+    public static void await(CountDownLatch latch) {
+        try {
+            latch.await();
+        }
+        catch (InterruptedException e) {
+            Utilities.exit("ERROR", e);
+        }
+    }
+
+    /**
+     * Shuts down and awaits termination of the specified
+     * {@code ExecutorService}. The Java Virtual Machine is terminated if an
+     * {@code InterruptedException} is thrown while awaiting termination
+     * of the executor service.
+     * @param es the executor service to be shut down
+     * @throws NullPointerException if {@code es == null}
+     */
+    public static void shutdownExecService(ExecutorService es) {
+        try {
+            es.shutdown();
+            es.awaitTermination(Long.MAX_VALUE, TimeUnit.DAYS);
+        }
+        catch (InterruptedException e) {
+            Utilities.exit("ERROR", e);
+        }
+    }
+
+}
+

blbutil/SampleFileIt.java

@@ -18,7 +18,7 @@
  */
 package blbutil;
 
-import beagleutil.Samples;
+import vcf.Samples;
 
 /**
  * <p>An iterator for records in a file.  Each records contains

blbutil/Utilities.java

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

bref/AsIsBref3Writer.java

@@ -19,11 +19,10 @@
 package bref;
 
 import beagleutil.ChromIds;
-import beagleutil.Samples;
 import blbutil.FileUtil;
+import blbutil.Utilities;
 import ints.IntArray;
 import ints.IntList;
-import blbutil.Utilities;
 import java.io.ByteArrayOutputStream;
 import java.io.DataOutputStream;
 import java.io.File;
@@ -38,6 +37,7 @@ import java.util.List;
 import java.util.Set;
 import vcf.Marker;
 import vcf.RefGTRec;
+import vcf.Samples;
 
 /**
  * <p>Class {@code AsIsBref3Writer} writes VCF data with phased, non-missing

bref/Bref3.java

@@ -18,7 +18,6 @@
  */
 package bref;
 
-import beagleutil.Samples;
 import blbutil.Const;
 import blbutil.FileIt;
 import blbutil.InputIt;
@@ -27,6 +26,7 @@ import blbutil.Utilities;
 import java.io.File;
 import vcf.RefGTRec;
 import vcf.RefIt;
+import vcf.Samples;
 
 /**
  * <p>Class {@code Bref3} converts files in VCF format into 

bref/Bref3It.java

@@ -18,7 +18,6 @@
  */
 package bref;
 
-import beagleutil.Samples;
 import blbutil.FileUtil;
 import blbutil.Filter;
 import blbutil.SampleFileIt;
@@ -33,6 +32,7 @@ import java.util.Deque;
 import java.util.NoSuchElementException;
 import vcf.Marker;
 import vcf.RefGTRec;
+import vcf.Samples;
 
 /**
  * <p>Class {@code Bref3It} represents  an iterator whose {@code next()} which

bref/Bref3Reader.java

@@ -19,13 +19,12 @@
 package bref;
 
 import beagleutil.ChromIds;
-import beagleutil.Samples;
-import ints.CharArray;
 import blbutil.Const;
 import blbutil.Filter;
+import blbutil.Utilities;
+import ints.CharArray;
 import ints.IntArray;
 import ints.UnsignedByteArray;
-import blbutil.Utilities;
 import java.io.DataInput;
 import java.io.IOException;
 import java.util.ArrayList;
@@ -35,6 +34,7 @@ import java.util.List;
 import vcf.BasicMarker;
 import vcf.Marker;
 import vcf.RefGTRec;
+import vcf.Samples;
 import vcf.SeqCodedRefGTRec;
 
 /**
@@ -82,8 +82,10 @@ public final class Bref3Reader {
         } catch (IOException ex) {
             Utilities.exit(READ_ERR, ex);
         }
+        boolean[] isDiploid = new boolean[sampleIds.length];
+        Arrays.fill(isDiploid, true);
         this.program = programString;
-        this.samples = Samples.fromIds(sampleIds);
+        this.samples = Samples.fromIds(sampleIds, isDiploid);
         this.markerFilter = markerFilter;
         this.nHaps = 2*samples.nSamples();
         this.byteBuffer = new byte[2*nHaps];

bref/BrefWriter.java

@@ -18,9 +18,9 @@
  */
 package bref;
 
-import beagleutil.Samples;
 import java.io.Closeable;
 import vcf.RefGTRec;
+import vcf.Samples;
 
 /**
  * <p>Interface {@code BrefWrites} writes phased, non-missing genotypes to a

bref/CompressBref3Writer.java

@@ -18,11 +18,11 @@
  */
 package bref;
 
-import beagleutil.Samples;
 import java.io.File;
 import java.util.ArrayList;
 import java.util.List;
 import vcf.RefGTRec;
+import vcf.Samples;
 
 /**
  * <p>Class {@code CompressBref3Writer} writes phased, non-missing genotypes

bref/SeqCoder3.java

@@ -18,7 +18,6 @@
  */
 package bref;
 
-import beagleutil.Samples;
 import ints.IntArray;
 import ints.IntList;
 import java.util.ArrayList;
@@ -27,11 +26,12 @@ import java.util.Collections;
 import java.util.List;
 import vcf.Marker;
 import vcf.RefGTRec;
+import vcf.Samples;
 import vcf.SeqCodedRefGTRec;
 
 /**
  * <p>Class {@code SeqCoder3} compresses a sequence of allele-coded
- * {@code RefGTRec} objects.  The class is designed for use with brev v3 format.
+ * {@code RefGTRec} objects.  The class is designed for use with bref v3 format.
  * Compression is performed by storing the list of distinct allele sequences
  * and the allele sequence carried by each haplotype.
  * </p>

debian/changelog

+beagle (5.1-191108+dfsg-1) UNRELEASED; urgency=medium
+
+  * New upstream release.
+
+ -- Dylan Aïssi <daissi@debian.org>  Wed, 20 Nov 2019 08:00:44 +0100
+
 beagle (5.1-190921+dfsg-2) unstable; urgency=medium
 
   * Fix watch file.

imp/HaplotypeCoder.java

@@ -18,8 +18,8 @@
  */
 package imp;
 
-import ints.IntArray;
 import ints.IndexArray;
+import ints.IntArray;
 import java.util.Arrays;
 import vcf.GT;
 import vcf. RefGT;

imp/ImpData.java

@@ -18,20 +18,20 @@
  */
 package imp;
 
-import beagleutil.Samples;
+import ints.IndexArray;
 import ints.IntArray;
 import ints.IntList;
-import ints.IndexArray;
 import java.util.Arrays;
 import java.util.stream.IntStream;
-import vcf.GeneticMap;
 import main.Par;
 import vcf.Data;
 import vcf.GT;
+import vcf.GeneticMap;
 import vcf.MarkerIndices;
 import vcf.Markers;
 import vcf.RefGT;
 import vcf. RefGTRec;
+import vcf.Samples;
 
 /**
  * <p>Class {@code ImpData} contains the input data for imputation of
@@ -61,6 +61,7 @@ public class ImpData {
     private final int nClusters;
     private final int nRefHaps;
     private final int nTargHaps;
+    private final int nInputTargHaps;
     private final int nHaps;
 
     /**
@@ -100,6 +101,11 @@ public class ImpData {
         this.nClusters = targClustStartEnd.length - 1;
         this.nRefHaps = refGT.nHaps();
         this.nTargHaps = phasedTarg.nHaps();
+        Samples targSamples = phasedTarg.samples();
+        this.nInputTargHaps = IntStream.range(0, targSamples.nSamples())
+                .parallel()
+                .map(j -> (targSamples.isDiploid(j) ? 2 : 1))
+                .sum();
         this.nHaps = refGT.nHaps() + phasedTarg.nHaps();
         this.hapToSeq = hapToSeq(data.restrictRefGT(), phasedTarg, targClustStartEnd);
         this.refClusterStart = refClustStart(targClustStartEnd, targToRef);
@@ -370,6 +376,15 @@ public class ImpData {
         return nTargHaps;
     }
 
+    /**
+     * Return the number of input target haplotypes. Each sample with
+     * haploid input data contributes only one haplotype to the total.
+     * @return the number of input target haplotypes
+     */
+    public int nInputTargHaps() {
+        return nInputTargHaps;
+    }
+
     /**
      * Returns the specified target marker cluter alleles for the
      * reference and target haplotypes.  Alleles for the

imp/ImpIbs.java

@@ -18,10 +18,10 @@
  */
 package imp;
 
-import ints.IntArray;
-import ints.IntList;
 import blbutil.Utilities;
 import ints.IndexArray;
+import ints.IntArray;
+import ints.IntList;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.List;

imp/ImpStates.java

@@ -18,12 +18,12 @@
  */
 package imp;
 
-import main.CompHapSegment;
 import ints.IntIntMap;
 import ints.IntList;
 import java.util.PriorityQueue;
 import java.util.Random;
 import java.util.stream.IntStream;
+import main.CompHapSegment;
 
 /**
  * <p>Class {@code ImpStates} identifies a list of pseudo-reference haplotypes
@@ -165,23 +165,22 @@ public final class ImpStates {
     }
 
     private int copyData(int targHap, int[][] hapIndices, boolean[][] alMatch) {
-        int nCompositeHaps = q.size();
-        initializeCopy(nCompositeHaps);
+        int nCompHaps = q.size();
         int shiftedTargHap = impData.nRefHaps() + targHap;
+        initializeCopy(nCompHaps);
         for (int m=0; m<nClusters; ++m) {
-            int targAl = impData.allele(m, shiftedTargHap);
-            for (int j=0; j<nCompositeHaps; ++j) {
+            int targAllele = impData.allele(m, shiftedTargHap);
+            for (int j=0; j<nCompHaps; ++j) {
                 if (m==compHapToEnd[j]) {
                     ++compHapToListIndex[j];
                     compHapToHap[j] = compositeHapToHap[j].get(compHapToListIndex[j]);
                     compHapToEnd[j] = compositeHapToEnd[j].get(compHapToListIndex[j]);
-                    assert compHapToHap[j] < impData.nRefHaps();
                 }
                 hapIndices[m][j] = compHapToHap[j];
-                alMatch[m][j] = impData.allele(m, compHapToHap[j])==targAl;
+                alMatch[m][j] = impData.allele(m, compHapToHap[j])==targAllele;
             }
         }
-        return nCompositeHaps;
+        return nCompHaps;
     }
 
     private void initializeCopy(int nSlots) {

imp/ImputedRecBuilder.java

@@ -53,7 +53,7 @@ public final class ImputedRecBuilder {
 
     private final Marker marker;
     private final int nAlleles;
-    private final int nHaps;
+    private final int nInputTargHaps;
     private final boolean ap;
     private final boolean gp;
     private final float[] sumAlProbs;
@@ -61,35 +61,36 @@ public final class ImputedRecBuilder {
     private final String[] homRefField;
     private final StringBuilder sampleData;
 
-    private int sampleCnt;
+    private int hapCnt;
 
     /**
      * Constructs a new {@code ImputedRecBuilder} instance for the specified
      * number of samples.
      *
      * @param marker the marker corresponding to the VCF record
-     * @param nSamples the number of samples
+     * @param nInputTargHaps the number of input target haplotypes for haploid
+     * and diploid samples
      * @param ap {@code true} if posterior allele probabilities are to be printed
      * @param gp {@code true} if posterior genotype probabilities are to be printed
-     * @throws IllegalArgumentException if {@code nSamples < 1}
+     * @throws IllegalArgumentException if {@code nInputTargHaps < 1}
      * @throws NullPointerException if {@code marker == null}
      */
-    public ImputedRecBuilder(Marker marker, int nSamples, boolean ap,
+    public ImputedRecBuilder(Marker marker, int nInputTargHaps, boolean ap,
             boolean gp) {
-        if (nSamples < 1) {
-            throw new IllegalArgumentException(String.valueOf(nSamples));
+        if (nInputTargHaps < 1) {
+            throw new IllegalArgumentException(String.valueOf(nInputTargHaps));
         }
         this.marker = marker;
         this.nAlleles = marker.nAlleles();
-        this.nHaps = 2*nSamples;
+        this.nInputTargHaps = nInputTargHaps;
         this.ap = ap;
         this.gp = gp;
         this.sumAlProbs = new float[nAlleles];
         this.sumAlProbs2 = new float[nAlleles];
-        this.sampleData = new StringBuilder(200 + nSamples*10);
+        this.sampleData = new StringBuilder(200 + nInputTargHaps*5);
         this.homRefField = (ap || gp) ? homRefFields(ap, gp)
                 : DEFAULT_HOM_REF_FIELDS;
-        this.sampleCnt = 0;
+        this.hapCnt = 0;
     }
 
     /**
@@ -101,21 +102,23 @@ public final class ImputedRecBuilder {
     }
 
     /**
-     * Returns the required number of samples in the VCF record.
-     * @return the required number of samples in the VCF record
+     * Returns the number of input target haplotypes for haploid
+     * and diploid samples.
+     * @return the number of input target haplotypes for haploid
+     * and diploid samples
      */
-    public int nSamples() {
-        return nHaps/2;
+    public int nInputTargHaps() {
+        return nInputTargHaps;
     }
 
     /**
-     * Returns the number of times the {@code this.addSampleData()}
-     * method has been invoked.
-     * @return the number of times the {@code this.addSampleData()}
-     * method has been invoked
+     * Returns the number of imputed alleles added by the addSampleData()
+     * methods.
+     * @return the number of imputed alleles added by the addSampleData()
+     * methods
      */
-    public int sampleCnt() {
-        return sampleCnt;
+    public int hapCnt() {
+        return hapCnt;
     }
 
     /**
@@ -132,7 +135,7 @@ public final class ImputedRecBuilder {
      * @throws NullPointerException if {@code a1 == null || a2 == null}
      */
     public void addSampleData(float[] a1, float[] a2) {
-        ++sampleCnt;
+        hapCnt+=2;
         if (a1[0]==1.0f && a2[0]==1.0f && a1.length < DEFAULT_HOM_REF_FIELDS.length) {
             sampleData.append(homRefField[a1.length]);
         }
@@ -176,6 +179,39 @@ public final class ImputedRecBuilder {
         }
     }
 
+    /**
+     * Scales the specified probabilities for each allele to each sum to 1.0,
+     * and adds the sample data to the VCF record.  The contract
+     * for this method is undefined if any element of the specified arrays is
+     * not a finite non-negative number.
+     * @param a1 the allele probabilities
+     * @throws IndexOutOfBoundsException if
+     * {@code a1.length < this.marker().nAlleles()}
+     * @throws IndexOutOfBoundsException if
+     * {@code a2.length < this.marker().nAlleles()}
+     * @throws NullPointerException if {@code a1 == null || a2 == null}
+     */
+    public void addSampleData(float[] a1) {
+        ++hapCnt;
+        scale(a1);
+        sampleData.append(Const.tab);
+        sampleData.append(maxIndex(a1));
+        for (int a=1; a<nAlleles; ++a) {
+            float dose = a1[a];
+            float dose2 = a1[a]*a1[a];
+            sumAlProbs[a] += dose;
+            sumAlProbs2[a] += dose2;
+            sampleData.append( (a==1) ? Const.colon : Const.comma );
+            sampleData.append(DS_VALS[(int) Math.rint(100*dose)]);
+        }
+        if (ap) {
+            for (int a=1; a<nAlleles; ++a) {
+                sampleData.append( (a==1) ? Const.colon : Const.comma );
+                sampleData.append(DS_VALS[(int) Math.rint(100*a1[a])]);
+            }
+        }
+    }
+
     private static void scale(float[] fa) {
         float sum = 0f;
         for (float f : fa) {
@@ -205,11 +241,11 @@ public final class ImputedRecBuilder {
       *@param out the {@code PrintWriter} to which the VCF record will be
      * printed
      * @throws IllegalStateException if
-     * {@code this.nSamples() != this.sampleCnt()}
+     * {@code this.hapCnt() != this.nInputTargHaps()}
      * @throws NullPointerException if {@code out == null}
      */
     public void printRec(PrintWriter out, boolean isImputed) {
-        if (2*sampleCnt != nHaps) {
+        if (hapCnt != nInputTargHaps) {
             throw new IllegalStateException("inconsistent data");
         }
         printMarkerFields(marker, out);
@@ -273,7 +309,7 @@ public final class ImputedRecBuilder {
             }
             for (int a=1; a<nAlleles; ++a) {
                 out.print( (a==1) ? ";AF=" : Const.comma);
-                out.print(DF4.format(sumAlProbs[a]/nHaps));
+                out.print(DF4.format(sumAlProbs[a]/nInputTargHaps));
             }
             if (marker.end()!=-1) {
                 out.print(";END=");
@@ -292,7 +328,7 @@ public final class ImputedRecBuilder {
         }
         else {
             float sum2 = sumAlProbs2[allele];
-            float meanTerm = sum*sum/(nHaps);
+            float meanTerm = sum*sum/(nInputTargHaps);
             float num = (sum2 - meanTerm);
             float den = (sum - meanTerm);
             return num <= 0 ? 0f : num/den;

imp/ImputedVcfWriter.java

@@ -24,9 +24,10 @@ import java.io.PrintWriter;
 import java.util.Arrays;
 import java.util.concurrent.atomic.AtomicReferenceArray;
 import java.util.stream.IntStream;
+import vcf.MarkerIndices;
 import vcf.Markers;
 import vcf. RefGT;
-import vcf.MarkerIndices;
+import vcf.Samples;
 
 /**
  * <p>Class {@code ImputedVcfWriter} writes observed and imputed genotypes
@@ -40,6 +41,7 @@ import vcf.MarkerIndices;
 public class ImputedVcfWriter {
 
     private final ImpData impData;
+    private final Samples targSamples;
     private final int targCluster;
     private final int refStart;
     private final int clustEnd;
@@ -74,6 +76,7 @@ public class ImputedVcfWriter {
             throw new IndexOutOfBoundsException(String.valueOf(refEnd));
         }
         this.impData = impData;
+        this.targSamples = impData.targSamples();
         this.targCluster = targCluster;
         if (targCluster==0) {
             this.refStart = refStart;
@@ -108,15 +111,16 @@ public class ImputedVcfWriter {
         }
         RefHapHash refHapHash = new RefHapHash(stateProbs, targCluster,
                 impData.refGT(), refStart, refEnd);
-
         ImputedRecBuilder[] recBuilders = recBuilders();
         float[][] a1Probs = alProbs();
         float[][] a2Probs = alProbs();
         boolean[] isImputed = isImputed();
         for (int h=0, n=stateProbs.length(); h<n; h+=2) {
+            boolean isDiploid = targSamples.isDiploid(h>>1);
             setAlProbs(stateProbs.get(h), refHapHash, a1Probs);
             setAlProbs(stateProbs.get(h+1), refHapHash, a2Probs);
 
+            if (isDiploid) {
                 for (int m=0; m<a1Probs.length; ++m) {
                     if (isImputed[m]==false) {
                         setToObsAlleles(a1Probs, a2Probs, m, h);
@@ -126,6 +130,17 @@ public class ImputedVcfWriter {
                     Arrays.fill(a2Probs[m], 0f);
                 }
             }
+            else {
+                for (int m=0; m<a1Probs.length; ++m) {
+                    if (isImputed[m]==false) {
+                        setToObsAlleles(a1Probs, a2Probs, m, h);
+                    }
+                    recBuilders[m].addSampleData(a1Probs[m]);
+                    Arrays.fill(a1Probs[m], 0f);
+                    Arrays.fill(a2Probs[m], 0f);
+                }
+            }
+        }
         for (int m=0; m<a1Probs.length; ++m) {
             recBuilders[m].printRec(out, isImputed[m]);
         }
@@ -203,10 +218,10 @@ public class ImputedVcfWriter {
         RefGT refGT = impData.refGT();
         boolean gp = impData.par().gp();
         boolean ap = impData.par().ap();
-        int nTargSamples = impData.nTargSamples();
+        int nInputTargHaps = impData.nInputTargHaps();
         return IntStream.range(refStart, refEnd)
                 .mapToObj(m -> new ImputedRecBuilder(refGT.marker(m),
-                        nTargSamples, ap, gp))
+                        nInputTargHaps, ap, gp))
                 .toArray(ImputedRecBuilder[]::new);
     }
 

main/Main.java

@@ -18,7 +18,6 @@
  */
 package main;
 
-import vcf.GeneticMap;
 import beagleutil.ChromInterval;
 import blbutil.Const;
 import blbutil.FileIt;
@@ -26,31 +25,32 @@ import blbutil.Filter;
 import blbutil.InputIt;
 import blbutil.SampleFileIt;
 import blbutil.Utilities;
+import bref.Bref3It;
 import imp.ImpData;
 import imp.ImpLS;
+import imp.StateProbs;
+import ints.LongArray;
 import java.io.File;
 import java.util.HashSet;
 import java.util.Locale;
 import java.util.Set;
-import vcf.AllData;
-import bref.Bref3It;
-import imp.StateProbs;
-import ints.LongArray;
 import java.util.concurrent.atomic.AtomicReferenceArray;
 import java.util.function.Supplier;
 import phase.FixedPhaseData;
-import vcf.VcfIt;
+import vcf.AllData;
 import vcf.Data;
-import vcf.IntervalVcfIt;
-import vcf.Marker;
 import vcf.FilterUtil;
 import vcf.GT;
-import vcf.Markers;
-import vcf.TargetData;
-import vcf.RefIt;
 import vcf.GTRec;
+import vcf.GeneticMap;
 import vcf.HapsGT;
+import vcf.IntervalVcfIt;
+import vcf.Marker;
+import vcf.Markers;
 import vcf.RefGTRec;
+import vcf.RefIt;
+import vcf.TargetData;
+import vcf.VcfIt;
 
 /**
  * Class {@code Main} is the entry class for the Beagle program.
@@ -65,20 +65,20 @@ public class Main {
      * The program name and version.
      */
     public static final String VERSION = "(version 5.1)";
-    public static final String PROGRAM = "beagle.21Sep19.ec3.jar";
-    public static final String COMMAND = "java -jar beagle.21Sep19.ec3.jar";
+    public static final String PROGRAM = "beagle.08Nov19.3ec.jar";
+    public static final String COMMAND = "java -jar beagle.08Nov19.3ec.jar";
 
     /**
      * The copyright string.
      */
-    public static final String COPYRIGHT = "Copyright (C) 2014-2019 Brian L. Browning";
+    public static final String COPYRIGHT = "Copyright (C) 2014-2018 Brian L. Browning";
 
     /**
      * The program name and a brief help message.
      */
     public static final String SHORT_HELP = Main.PROGRAM + " " + VERSION
             + Const.nl + Main.COPYRIGHT
-            + Const.nl + "Enter \"java -jar beagle.21Sep19.ec3.jar\" to "
+            + Const.nl + "Enter \"java -jar beagle.08Nov19.3ec.jar\" to "
             + "list command line argument";
 
     private final Par par;
@@ -108,8 +108,7 @@ public class Main {
         runStats.printStartInfo();
 
         try (Data data = data(par, runStats);
-                WindowWriter winOut = new WindowWriter(data.targGT().samples(),
-                        par.out())) {
+                WindowWriter winOut = new WindowWriter(par, data.targGT().samples())) {
             Main main = new Main(par, data, winOut, runStats);
             main.phaseData();
             runStats.printSummaryAndClose(data.nTargMarkersSoFar(),
@@ -255,8 +254,7 @@ public class Main {
                             + Const.nl);
                 }
                 FileIt<String> it = InputIt.fromGzipFile(par.ref());
-                refIt = RefIt.create(it, sFilter, mFilter2,
-                        RefIt.MAX_EM_BUFFER_SIZE);
+                refIt = RefIt.create(it, sFilter, mFilter2);
             }
             if (chromInt!=null) {
                 refIt = new IntervalVcfIt<>(refIt, chromInt);

main/MainHelper.java

@@ -18,11 +18,11 @@
  */
 package main;
 
-import phase.EstPhase;
-import phase.PhaseLS;
 import java.util.Random;
+import phase.EstPhase;
 import phase.FixedPhaseData;
 import phase.PhaseData;
+import phase.PhaseLS;
 import vcf.GT;
 import vcf. RefGT;
 
@@ -107,7 +107,8 @@ public class MainHelper {
             int nIts = par.burnin() + par.iterations();
             PhaseData phaseData = new PhaseData(fpd, estPhase, recombFactor,
                     nIts, rand.nextLong());
-            return PhaseLS.runStage2(phaseData);
+            GT result = PhaseLS.runStage2(phaseData);
+            return result;
         }
         else {
             return estPhase.hapsGT();