[med-svn] [Git][med-team/beagle][upstream] New upstream version 220208
Andreas Tille (@tille)
gitlab at salsa.debian.org
Thu Feb 17 17:55:57 GMT 2022
Andreas Tille pushed to branch upstream at Debian Med / beagle
Commits:
3b5c4164 by Andreas Tille at 2022-02-17T17:19:37+01:00
New upstream version 220208
- - - - -
9 changed files:
- bref/Bref3.java
- bref/UnBref3.java
- main/Main.java
- main/Par.java
- main/RunStats.java
- phase/HmmParamData.java
- phase/ParamEstimates.java
- phase/PhaseData.java
- phase/PhaseLS.java
Changes:
=====================================
bref/Bref3.java
=====================================
@@ -38,7 +38,7 @@ import vcf.Samples;
*/
public class Bref3 {
- private static final String program = "bref3.28Jun21.220.jar";
+ private static final String program = "bref3.08Feb22.fa4.jar";
/**
* The {@code main()} method is the entry point to the bref program.
=====================================
bref/UnBref3.java
=====================================
@@ -37,7 +37,7 @@ import vcf.VcfWriter;
*/
public class UnBref3 {
- private static final String program = "unbref3.28Jun21.220.jar";
+ private static final String program = "unbref3.08Feb22.fa4.jar";
/**
* The {@code main()} method is the entry point to the bref program.
=====================================
main/Main.java
=====================================
@@ -68,17 +68,17 @@ public class Main {
/**
* The program name and version.
*/
- private static final String VERSION = "(version 5.2)";
+ private static final String VERSION = "(version 5.3)";
/**
- * The program name and version.
+ * The program name and commit version.
*/
- public static final String PROGRAM = "beagle.28Jun21.220.jar";
+ public static final String PROGRAM = "beagle.08Feb22.fa4.jar";
/**
* The command to invoke the program.
*/
- public static final String COMMAND = "java -jar beagle.28Jun21.220.jar";
+ public static final String COMMAND = "java -jar beagle.08Feb22.fa4.jar";
/**
* The copyright string.
@@ -191,8 +191,8 @@ public class Main {
while (pd.it()<nIts) {
long t0 = System.nanoTime();
PhaseLS.runStage1(pd);
- if (pd.it()==par.burnin()) {
- runStats.printEstimatedNe(pd.ne());
+ if (pd.it()==par.burnin() && par.em()) {
+ runStats.printEstimatedParameters(pd.ne(), pd.pMismatch());
}
runStats.printStage1Info(pd.it(), (System.nanoTime() - t0));
pd.incrementIt();
=====================================
main/Par.java
=====================================
@@ -35,6 +35,7 @@ import java.util.Map;
public final class Par {
private final String[] args;
+ private final boolean noNThreads;
// data parameters
private final File gt;
@@ -64,6 +65,7 @@ public final class Par {
private final boolean gp;
// general parameters
+ private final boolean em;
private final float ne;
private final float err;
private final float window;
@@ -93,7 +95,8 @@ public final class Par {
private static final boolean D_GP = false;
// default general parameters
- private static final int D_NE = 1_000_000;
+ private static final boolean D_EM = true;
+ private static final int D_NE = 100_000;
private static final float D_ERR = -Float.MIN_VALUE;
private static final float D_WINDOW = 40.0f;
private static final float D_OVERLAP = 2.0f;
@@ -150,17 +153,20 @@ public final class Par {
imp_step = Validate.floatArg("imp-step", argsMap, false, D_IMP_STEP, FMIN, FMAX);
imp_nsteps = Validate.intArg("imp-nsteps", argsMap, false, D_IMP_NSTEPS, 1, IMAX);
cluster = Validate.floatArg("cluster", argsMap, false, D_CLUSTER, 0.0f, FMAX);
- ap = Validate.booleanArg("ap", argsMap, D_AP, false);
- gp = Validate.booleanArg("gp", argsMap, D_GP, false);
+ ap = Validate.booleanArg("ap", argsMap, false, D_AP);
+ gp = Validate.booleanArg("gp", argsMap, false, D_GP);
// general parameters
+ em = Validate.booleanArg("em", argsMap, false, D_EM);
ne = Validate.floatArg("ne", argsMap, false, D_NE, FMIN, FMAX);
err = Validate.floatArg("err", argsMap, false, D_ERR, -FMIN, FMAX);
window = Validate.floatArg("window", argsMap, false, D_WINDOW, FMIN, 1000);
overlap = Validate.floatArg("overlap", argsMap, false, D_OVERLAP, FMIN, IMAX);
buffer = Validate.floatArg("buffer", argsMap, false, D_BUFFER, FMIN, IMAX);
seed = Validate.longArg("seed", argsMap, false, D_SEED, LMIN, LMAX);
- nthreads = modNthreads(Validate.intArg("nthreads", argsMap, false, D_NTHREADS, 1, IMAX));
+ int rawNThreads = Validate.intArg("nthreads", argsMap, false, D_NTHREADS, 1, IMAX);
+ noNThreads = rawNThreads==D_NTHREADS;
+ nthreads = noNThreads ? Runtime.getRuntime().availableProcessors() : rawNThreads;
// undocumented parameters
truth = Validate.getFile(Validate.stringArg("truth", argsMap, false, null, null));
@@ -176,6 +182,16 @@ public final class Par {
return args.clone();
}
+ /**
+ * Returns {@code true} if the command line does not include an
+ * nthreads parameter and returns {@code false otherwise}.
+ * @return {@code true} if the command line does not include an
+ * nthreads parameter
+ */
+ public boolean noNThreads() {
+ return noNThreads;
+ }
+
/**
* Returns a description of the Beagle command line arguments.
* @return a description of the Beagle command line arguments
@@ -209,6 +225,7 @@ public final class Par {
+ "general parameters ..." + nl
+ " ne=<effective population size> (default=" + D_NE + ")" + nl
+ " window=<window length in cM> (default=" + D_WINDOW + ")" + nl
+ + " em=<estimate ne and err parameters (true/false)> (default=" + D_EM + ")" + nl
+ " overlap=<window overlap in cM> (default=" + D_OVERLAP + ")" + nl
+ " seed=<random seed> (default=" + D_SEED + ")" + nl
+ " nthreads=<number of threads> (default: machine dependent)" + nl + nl;
@@ -223,21 +240,6 @@ public final class Par {
return ci;
}
- /**
- * Returns the nthreads parameter, which is equal to
- * {@code Runtime.getRuntime().availableProcessors()} if
- * {@code nthreads == Integer.MAX_VALUE}.
- * @return the nthreads parameter
- */
- private static int modNthreads(int nthreads) {
- if (nthreads==Integer.MAX_VALUE) {
- return Runtime.getRuntime().availableProcessors();
- }
- else {
- return nthreads;
- }
- }
-
// data parameters
/**
@@ -430,6 +432,14 @@ public final class Par {
// general parameters
+ /**
+ * Returns the em parameter.
+ * @return the em parameter
+ */
+ public boolean em() {
+ return em;
+ }
+
/**
* Returns the ne parameter
* @return the ne parameter
=====================================
main/RunStats.java
=====================================
@@ -23,6 +23,7 @@ import blbutil.FileUtil;
import blbutil.Utilities;
import java.io.File;
import java.io.PrintWriter;
+import java.util.Arrays;
import java.util.Optional;
import vcf.Data;
import vcf.Marker;
@@ -62,9 +63,16 @@ public class RunStats {
* file and to standard output.
*/
public void printStartInfo() {
+ String[] argList = par.args();
+ if (par.noNThreads()) {
+ // add nthreads parameter to list of command line parameters
+ int oldLength = argList.length;
+ argList = Arrays.copyOf(argList, argList.length+1);
+ argList[oldLength] = "nthreads=" + par.nthreads();
+ }
Utilities.duoPrint(log, Main.SHORT_HELP + Const.nl);
Utilities.duoPrintln(log, "Start time: " + Utilities.timeStamp());
- Utilities.duoPrint(log, Utilities.commandLine(Main.PROGRAM, par.args()));
+ Utilities.duoPrint(log, Utilities.commandLine(Main.PROGRAM, argList));
if (par.ped() != null) {
String s = Const.nl + "WARNING: This version will not model"
+ " duos or trios in the pedigree file";
@@ -253,15 +261,14 @@ public class RunStats {
/**
* Prints the specified estimated effective population size.
* @param ne the estimated effective population size
- * @throws IllegalArgumentException if {@code ne < 1}
+ * @param pMismatch the estimated allele mismatch parameter
*/
- public void printEstimatedNe(long ne) {
- if (ne<1) {
- throw new IllegalArgumentException(String.valueOf(ne));
- }
+ public void printEstimatedParameters(long ne, float pMismatch) {
Utilities.duoPrintln(log, "");
Utilities.duoPrint(log, String.format("%1$-31s", "Estimated ne:"));
Utilities.duoPrintln(log, String.valueOf(ne));
+ Utilities.duoPrint(log, String.format("%1$-31s", "Estimated err:"));
+ Utilities.duoPrintln(log, String.format("%1$7.1e", pMismatch));
}
/**
=====================================
phase/HmmParamData.java
=====================================
@@ -46,9 +46,6 @@ public class HmmParamData {
private final float[][] savedBwd;
private final float[] emProbs;
- private final boolean[] isHom;
- private final float maxErrMaf;
- private final FloatArray markerToMaf;
private int mismatchCnt;
private double sumMismatchProb;
private double sumGenDist;
@@ -58,18 +55,9 @@ public class HmmParamData {
* Creates a {@code HmmParamData} instance for the specified data.
*
* @param phaseIbs the IBS haplotype segments
- * @param maxErrMaf the maximum minor allele frequency for which allele
- * allele mismatch probabilities will be estimated
- * @throws IllegalArgumentException if
- * {@code maxErrMaf <= 0.0 || maxErrMaf >= 1.0
- * || Float.isFinite(maxErrMaf) == false}
* @throws NullPointerException if {@code phaseIbs == null}
*/
- public HmmParamData(PbwtPhaseIbs phaseIbs, float maxErrMaf) {
- if (maxErrMaf<=0.0 || maxErrMaf>=1.0
- || Float.isFinite(maxErrMaf)==false) {
- throw new IllegalArgumentException(String.valueOf(maxErrMaf));
- }
+ public HmmParamData(PbwtPhaseIbs phaseIbs) {
this.phaseData = phaseIbs.phaseData();
FixedPhaseData fpd = phaseData.fpd();
this.gt = fpd.stage1TargGT();
@@ -86,9 +74,6 @@ public class HmmParamData {
float pMismatch = phaseData.pMismatch();
this.emProbs = new float[] {1f - pMismatch, pMismatch};
- this.isHom = new boolean[nMarkers];
- this.markerToMaf = fpd.stage1Maf();
- this.maxErrMaf = maxErrMaf;
this.mismatchCnt = 0;
this.sumMismatchProb = 0.0;
this.sumGenDist = 0.0;
@@ -113,6 +98,10 @@ public class HmmParamData {
public void addEstimationData(ParamEstimates paramEst) {
paramEst.addMismatchData(mismatchCnt, sumMismatchProb);
paramEst.addSwitchData(sumGenDist, sumSwitchProb);
+ mismatchCnt = 0;
+ sumMismatchProb = 0.0;
+ sumGenDist = 0.0;
+ sumSwitchProb = 0.0;
}
/**
@@ -135,7 +124,6 @@ public class HmmParamData {
* {@code sample < 0 || sample >= this.nTargSamples()}
*/
public void update(int sample) {
- setIsHom(sample);
int nStates = states.ibsStates(sample, alMatch);
if (nStates>1) { // otherwise hFactor = nStates/(nStates - 1) is NaN
getParamData(alMatch[0], nStates);
@@ -183,10 +171,8 @@ public class HmmParamData {
mismatchSum += stateProb;
}
}
- if (isHom[m] && markerToMaf.get(m)<=maxErrMaf) {
- ++mismatchCnt;
- sumMismatchProb += (mismatchSum/stateSum);
- }
+ ++mismatchCnt;
+ sumMismatchProb += (mismatchSum/stateSum);
double switchProb = hFactor*(1.0f - jointStateSum/stateSum);
if (switchProb>0.0) {
sumGenDist += genDist.get(m);
@@ -194,12 +180,4 @@ public class HmmParamData {
}
return fwdSum;
}
-
- private void setIsHom(int sample) {
- for (int m=0; m<isHom.length; ++m) {
- int a1 = gt.allele1(m, sample);
- int a2 = gt.allele2(m, sample);
- isHom[m] = (a1>=0 && a1==a2);
- }
- }
}
=====================================
phase/ParamEstimates.java
=====================================
@@ -44,60 +44,55 @@ public class ParamEstimates {
}
/**
- * Records the specified allele mismatch data if {@code markerCnt > 0}
+ * Records the specified allele mismatch data if {@code markerCnt} and
+ * {@code pMismatchSum} are finite positive values.
* @param markerCnt the number of markers with mismatch data
* @param pMismatchSum the sum of estimated allele mismatch probabilities
*/
public void addMismatchData(int markerCnt, double pMismatchSum) {
- if (markerCnt>0) {
+ if (markerCnt>0 && pMismatchSum>0 && Double.isFinite(pMismatchSum)) {
mismatchData.add(new MismatchData(markerCnt, pMismatchSum));
}
}
/**
* Records the specified genetic distance and switch probability
- * if {@code genDistances > 0}.
+ * if {@code genDistances} and {@code switchProbs} are finite positive
+ * values.
* @param genDistances the list of genetic distance
* @param switchProbs the list of haplotype switch probabilities
*/
public void addSwitchData(double genDistances, double switchProbs) {
- if (genDistances>0) {
+ if (genDistances>0 && switchProbs>0
+ && Double.isFinite(genDistances) && Double.isFinite(switchProbs)) {
switchData.add(new RecombData(genDistances, switchProbs));
}
}
/**
- * Returns the estimated allele mismatch rate. Returns 0 if there is
- * no data to estimate the allele mismatch rate. The returned value is not
- * an atomic snapshot. Invocation in the absence of concurrent
+ * Returns the estimated allele mismatch rate. Returns {@code Float.NaN}
+ * if there is no data to estimate the allele mismatch rate. The returned
+ * value is not an atomic snapshot. Invocation in the absence of concurrent
* update will return an accurate result, but concurrent updates that
* occur white the sum is being calculated might not be incorporated
* in the returned result.
* @return the estimated allele mismatch rate
*/
public float pMismatch() {
- MismatchData[] mda = mismatchData.toArray(new MismatchData[0]);
- int sumMarkers = 0;
- double sumPMismatch = 0.0;
+ MismatchData[] mda = mismatchData.stream()
+ .sorted() // ensures sum of values is repeatable
+ .toArray(MismatchData[]::new);
+ long sumMarkers = 0;
+ double sumPMismatch = 0d;
for (MismatchData md : mda) {
sumMarkers += md.markerCnt;
sumPMismatch += md.pMismatchSum;
}
- return sumMarkers==0 ? 0f : (float) (sumPMismatch/sumMarkers);
+ return sumMarkers==0 ? Float.NaN : (float) (sumPMismatch/sumMarkers);
}
/**
- * Clears all data that has been added via the
- * {@code this.addMismatchData} and {@code this.addSwitchData()}
- * methods.
- */
- public void clear() {
- switchData.clear();
- mismatchData.clear();
- }
-
- /**
- * Returns the estimated recombination intensities. Returns {@code 0f}
+ * Returns the estimated recombination intensities. Returns {@code Float.NaN}
* if there is no data from which to estimate recombination intensities.
* The returned value is NOT an atomic snapshot. An accurate result is
* guaranteed only if no concurrent updates occur during method
@@ -105,17 +100,29 @@ public class ParamEstimates {
* @return the estimated recombination intensities
*/
public float recombIntensity() {
- RecombData[] rda = switchData.toArray(new RecombData[0]);
- double sumSwitches = 0;
- double sumDistances = 0;
+ RecombData[] rda = switchData.stream()
+ .sorted() // ensures sum of values is repeatable
+ .toArray(RecombData[]::new);
+ double sumSwitches = 0d;
+ double sumDistances = 0d;
for (RecombData rd : rda) {
sumSwitches += rd.switchProb;
sumDistances += rd.genDistance;
}
- return sumDistances==0d ? 0f : (float) (sumSwitches/sumDistances);
+ return sumDistances==0d ? Float.NaN : (float) (sumSwitches/sumDistances);
+ }
+
+ /**
+ * Clears all data that has been added via the
+ * {@code this.addMismatchData} and {@code this.addSwitchData()}
+ * methods.
+ */
+ public void clear() {
+ switchData.clear();
+ mismatchData.clear();
}
- private static class RecombData {
+ private static class RecombData implements Comparable<RecombData> {
private final double genDistance;
private final double switchProb;
@@ -124,9 +131,15 @@ public class ParamEstimates {
this.genDistance = genDistance;
this.switchProb = switchProb;
}
+
+ @Override
+ public int compareTo(RecombData o) {
+ int val = Double.compare(this.genDistance, o.genDistance);
+ return val!=0 ? val : Double.compare(this.switchProb, o.switchProb);
+ }
}
- private static class MismatchData {
+ private static class MismatchData implements Comparable<MismatchData> {
private final int markerCnt;
private final double pMismatchSum;
@@ -135,5 +148,11 @@ public class ParamEstimates {
this.markerCnt = markerCnt;
this.pMismatchSum = pMismatchSum;
}
+
+ @Override
+ public int compareTo(MismatchData o) {
+ int val = Double.compare(this.pMismatchSum, o.pMismatchSum);
+ return val!=0 ? val : Integer.compare(this.markerCnt, o.markerCnt);
+ }
}
}
=====================================
phase/PhaseData.java
=====================================
@@ -109,7 +109,7 @@ public class PhaseData {
* @return the effective population size
*/
private long ne(float recombIntensity) {
- return (long) Math.floor(25*recombIntensity*estPhase.fpd().nHaps());
+ return (long) Math.ceil(25*recombIntensity*estPhase.fpd().nHaps());
}
/**
=====================================
phase/PhaseLS.java
=====================================
@@ -18,7 +18,6 @@
*/
package phase;
-import blbutil.FloatArray;
import blbutil.MultiThreadUtils;
import blbutil.Utilities;
import java.util.Arrays;
@@ -38,8 +37,6 @@ import vcf.GT;
*/
public class PhaseLS {
- private static final Float MIN_PMISMATCH = 1e-12f;
-
private PhaseLS() {
// private constructor to prevent instantiation
}
@@ -51,16 +48,19 @@ public class PhaseLS {
*/
public static void runStage1(PhaseData pd) {
Par par = pd.fpd().par();
+ int it = pd.it();
int nThreads = par.nthreads();
int nSamples = pd.fpd().targGT().nSamples();
int nBurninIts = par.burnin();
PbwtPhaseIbs phaseIbs = pbwtPhaseIbs(pd);
- if (pd.it()==0) {
- initializeParameters(phaseIbs);
- }
- else if (pd.it()<nBurninIts) {
- resetPMismatchToDefault(pd); // only update pMismatch once
- updateParameters(phaseIbs);
+ if (par.em()) {
+ Random rand = new Random(pd.seed());
+ if (it==0) {
+ initializeParameters(phaseIbs, rand);
+ }
+ else if (it<nBurninIts) {
+ updateParameters(phaseIbs, rand);
+ }
}
AtomicInteger samples = new AtomicInteger(0);
ExecutorService es = Executors.newFixedThreadPool(nThreads);
@@ -87,100 +87,40 @@ public class PhaseLS {
return phaseIbs;
}
- private static void initializeParameters(PbwtPhaseIbs phaseIbs) {
+ private static void initializeParameters(PbwtPhaseIbs phaseIbs, Random rand) {
PhaseData pd = phaseIbs.phaseData();
- Random rand = new Random(pd.seed());
- ParamEstimates paramEst = new ParamEstimates();
- getParamEst(phaseIbs, paramEst, rand);
- if (isBigDecrease(pd, paramEst)) {
- decreaseRecombIntensity(phaseIbs, paramEst, rand);
- }
- else if (isBigIncrease(pd, paramEst)) {
- increaseRecombIntensity(phaseIbs, paramEst, rand);
- }
- int nTries = 0;
- float recombIntensity = paramEst.recombIntensity();
- boolean finished = recombIntensity==0f;
- while (++nTries<6 && finished==false) {
- pd.updateRecombIntensity(paramEst.recombIntensity());
- getParamEst(phaseIbs, paramEst, rand);
- recombIntensity = paramEst.recombIntensity();
- if (recombIntensity==0f) {
- finished = true;
- }
- else {
- boolean dec = paramEst.recombIntensity()<0.9f*pd.recombIntensity();
- boolean inc = paramEst.recombIntensity()>1.1111f*pd.recombIntensity();
- finished = dec==false && inc==false;
+ float prevRecInt = pd.recombIntensity();
+ int maxInitialIts = 15;
+ for (int j=0; j<maxInitialIts; ++j) {
+ updateParameters(phaseIbs, rand);
+ float recInt = pd.recombIntensity();
+ if (Math.abs(recInt - prevRecInt) <= 0.1*prevRecInt) {
+ break;
}
+ prevRecInt = recInt;
}
- updateParams(pd, paramEst);
}
- private static void updateParameters(PbwtPhaseIbs phaseIbs) {
+ private static void updateParameters(PbwtPhaseIbs phaseIbs, Random rand) {
+ ParamEstimates paramEst = getParamEst(phaseIbs, rand);
PhaseData pd = phaseIbs.phaseData();
- Random rand = new Random(phaseIbs.phaseData().seed());
- ParamEstimates paramEst = new ParamEstimates();
- getParamEst(phaseIbs, paramEst, rand);
- updateParams(pd, paramEst);
- }
-
- private static void updateParams(PhaseData pd, ParamEstimates paramEst) {
+ float prevPMismatch = pd.pMismatch();
float pMismatch = paramEst.pMismatch();
float recombIntensity = paramEst.recombIntensity();
- if (pMismatch>0f) {
- pd.updatePMismatch(Math.max(pMismatch, MIN_PMISMATCH));
+ if (Float.isFinite(pMismatch) && pMismatch>prevPMismatch) {
+ pd.updatePMismatch(pMismatch);
}
- if (recombIntensity>0f) {
+ if (Float.isFinite(recombIntensity) && recombIntensity>0f) {
pd.updateRecombIntensity(recombIntensity);
}
}
- private static void decreaseRecombIntensity(PbwtPhaseIbs phaseIbs,
- ParamEstimates paramEst, Random rand) {
- PhaseData pd = phaseIbs.phaseData();
- while (isBigDecrease(pd, paramEst)) {
- float nextRI = Math.min(0.333f*pd.recombIntensity(),
- paramEst.recombIntensity());
- pd.updateRecombIntensity(nextRI);
- getParamEst(phaseIbs, paramEst, rand);
- }
- }
-
- private static void increaseRecombIntensity(PbwtPhaseIbs phaseIbs,
- ParamEstimates paramEst, Random rand) {
- PhaseData pd = phaseIbs.phaseData();
- while (isBigIncrease(pd, paramEst)) {
- float nextRI = Math.max(3f*pd.recombIntensity(),
- paramEst.recombIntensity());
- pd.updateRecombIntensity(nextRI);
- getParamEst(phaseIbs, paramEst, rand);
- }
- }
-
- private static boolean isBigDecrease(PhaseData pd, ParamEstimates paramEst) {
- float estRI = paramEst.recombIntensity();
- return estRI>0f && estRI<0.8f*pd.recombIntensity();
- }
-
- private static boolean isBigIncrease(PhaseData pd, ParamEstimates paramEst) {
- float estRI = paramEst.recombIntensity();
- return estRI>0 && estRI>1.25f*pd.recombIntensity();
- }
-
- private static void resetPMismatchToDefault(PhaseData pd) {
- pd.updatePMismatch(Par.liStephensPMismatch(pd.fpd().nHaps()));
- }
-
- private static void getParamEst(PbwtPhaseIbs phaseIbs,
- ParamEstimates paramEst, Random rand) {
- paramEst.clear();
+ private static ParamEstimates getParamEst(PbwtPhaseIbs phaseIbs, Random rand) {
+ ParamEstimates paramEst = new ParamEstimates();
PhaseData pd = phaseIbs.phaseData();
FixedPhaseData fpd = pd.fpd();
int[] sampleIndices = samplesToAnalyze(pd, rand);
int nThreads = Math.min(fpd.par().nthreads(), sampleIndices.length);
- float mafPercentile = 0.05f;
- float maxErrMaf = maxMismatchMaf(fpd.stage1Maf(), mafPercentile);
double maxSum = 20000.0/nThreads;
int minIndices = 50;
AtomicInteger counter = new AtomicInteger(0);
@@ -188,14 +128,14 @@ public class PhaseLS {
for (int j=0; j<nThreads; ++j) {
es.submit(() -> {
try {
- HmmParamData hpd = new HmmParamData(phaseIbs, maxErrMaf);
+ HmmParamData hpd = new HmmParamData(phaseIbs);
int index = counter.getAndIncrement();
while ((hpd.sumSwitchProbs()<maxSum || index<minIndices)
&& index<sampleIndices.length) {
hpd.update(sampleIndices[index]);
index = counter.getAndIncrement();
+ hpd.addEstimationData(paramEst);
}
- hpd.addEstimationData(paramEst);
}
catch (Throwable t) {
Utilities.exit(t);
@@ -203,6 +143,7 @@ public class PhaseLS {
} );
}
MultiThreadUtils.shutdownExecService(es);
+ return paramEst;
}
private static int[] samplesToAnalyze(PhaseData pd, Random rand) {
@@ -220,23 +161,6 @@ public class PhaseLS {
}
}
- private static float maxMismatchMaf(FloatArray maf, float perc) {
- if (perc<=0f || perc>=1f || Float.isFinite(perc)==false) {
- throw new IllegalArgumentException(String.valueOf(perc));
- }
- double[] sortedMaf = IntStream.range(0, maf.size())
- .parallel()
- .mapToDouble(j -> maf.get(j))
- .sorted()
- .toArray();
- int start = Arrays.binarySearch(sortedMaf, Float.MIN_VALUE);
- if (start<0) {
- start = -start - 1;
- }
- int index = start + (int) Math.floor(perc * (sortedMaf.length - start));
- return (float) (index==sortedMaf.length ? Float.MIN_VALUE : sortedMaf[index]);
- }
-
/**
* Returns phased genotypes at all markers.
* @param pd estimated phased genotypes at stage 1 markers
View it on GitLab: https://salsa.debian.org/med-team/beagle/-/commit/3b5c416420dfba1f3b20402e52d8c90a5f13e89d
--
View it on GitLab: https://salsa.debian.org/med-team/beagle/-/commit/3b5c416420dfba1f3b20402e52d8c90a5f13e89d
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20220217/2ce3a37e/attachment-0001.htm>
More information about the debian-med-commit
mailing list