Commit 9ea91de8 authored by Dylan Aïssi's avatar Dylan Aïssi
Browse files

New upstream version 5.1-190824+dfsg

parent 7f4a872f
Loading
Loading
Loading
Loading
+3 −3
Original line number Diff line number Diff line
@@ -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);
    }
}
+7 −7
Original line number Diff line number Diff line
@@ -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;
        } ;
+190 −0
Original line number Diff line number Diff line
/*
 * 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();
            }
        }
    }
}
+92 −0
Original line number Diff line number Diff line
@@ -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);
    }
}
+13 −0
Original line number Diff line number Diff line
@@ -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
Loading