List of usage examples for org.apache.commons.math3.util MathArrays ebeSubtract
public static double[] ebeSubtract(double[] a, double[] b)
From source file:it.unibo.alchemist.model.implementations.positions.ContinuousGenericEuclidean.java
@Override public Position subtract(final Position other) { return new ContinuousGenericEuclidean(false, MathArrays.ebeSubtract(c, other.getCartesianCoordinates())); }
From source file:org.apache.solr.client.solrj.io.eval.EBESubtractEvaluator.java
@Override public Object doWork(Object first, Object second) throws IOException { if (null == first) { throw new IOException(String.format(Locale.ROOT, "Invalid expression %s - null found for the first value", toExpression(constructingFactory))); }/*from w w w.j a va 2s.co m*/ if (null == second) { throw new IOException(String.format(Locale.ROOT, "Invalid expression %s - null found for the second value", toExpression(constructingFactory))); } if (!(first instanceof List<?>)) { throw new IOException(String.format(Locale.ROOT, "Invalid expression %s - found type %s for the first value, expecting a list of numbers", toExpression(constructingFactory), first.getClass().getSimpleName())); } if (!(second instanceof List<?>)) { throw new IOException(String.format(Locale.ROOT, "Invalid expression %s - found type %s for the second value, expecting a list of numbers", toExpression(constructingFactory), first.getClass().getSimpleName())); } double[] result = MathArrays.ebeSubtract( ((List) first).stream().mapToDouble(value -> ((BigDecimal) value).doubleValue()).toArray(), ((List) second).stream().mapToDouble(value -> ((BigDecimal) value).doubleValue()).toArray()); List<Number> numbers = new ArrayList(); for (double d : result) { numbers.add(d); } return numbers; }
From source file:org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AlleleFrequencyCalculator.java
/** * Compute the probability of the alleles segregating given the genotype likelihoods of the samples in vc * * @param vc the VariantContext holding the alleles and sample information. The VariantContext * must have at least 1 alternative allele * @param refSnpIndelPseudocounts a total hack. A length-3 vector containing Dirichlet prior pseudocounts to * be given to ref, alt SNP, and alt indel alleles. Hack won't be necessary when we destroy the old AF calculators * @return result (for programming convenience) *//*from w w w .j a v a 2 s .c o m*/ @Override public AFCalculationResult getLog10PNonRef(final VariantContext vc, final int defaultPloidy, final int maximumAlternativeAlleles, final double[] refSnpIndelPseudocounts) { Utils.nonNull(vc, "vc is null"); final int numAlleles = vc.getNAlleles(); final List<Allele> alleles = vc.getAlleles(); Utils.validateArg(numAlleles > 1, "VariantContext has only a single reference allele, but getLog10PNonRef requires at least one at all " + vc); final double[] priorPseudocounts = alleles.stream().mapToDouble( a -> a.isReference() ? refPseudocount : (a.length() > 1 ? snpPseudocount : indelPseudocount)) .toArray(); double[] alleleCounts = new double[numAlleles]; final double flatLog10AlleleFrequency = -MathUtils.Log10Cache.get(numAlleles); // log10(1/numAlleles) double[] log10AlleleFrequencies = new IndexRange(0, numAlleles).mapToDouble(n -> flatLog10AlleleFrequency); double alleleCountsMaximumDifference = Double.POSITIVE_INFINITY; while (alleleCountsMaximumDifference > THRESHOLD_FOR_ALLELE_COUNT_CONVERGENCE) { final double[] newAlleleCounts = effectiveAlleleCounts(vc, log10AlleleFrequencies); alleleCountsMaximumDifference = Arrays.stream(MathArrays.ebeSubtract(alleleCounts, newAlleleCounts)) .map(Math::abs).max().getAsDouble(); alleleCounts = newAlleleCounts; final double[] posteriorPseudocounts = MathArrays.ebeAdd(priorPseudocounts, alleleCounts); // first iteration uses flat prior in order to avoid local minimum where the prior + no pseudocounts gives such a low // effective allele frequency that it overwhelms the genotype likelihood of a real variant // basically, we want a chance to get non-zero pseudocounts before using a prior that's biased against a variant log10AlleleFrequencies = new Dirichlet(posteriorPseudocounts).log10MeanWeights(); } double[] log10POfZeroCountsByAllele = new double[numAlleles]; double log10PNoVariant = 0; for (final Genotype g : vc.getGenotypes()) { if (!g.hasLikelihoods()) { continue; } final int ploidy = g.getPloidy() == 0 ? defaultPloidy : g.getPloidy(); final GenotypeLikelihoodCalculator glCalc = GL_CALCS.getInstance(ploidy, numAlleles); final double[] log10GenotypePosteriors = log10NormalizedGenotypePosteriors(g, glCalc, log10AlleleFrequencies); //the total probability log10PNoVariant += log10GenotypePosteriors[HOM_REF_GENOTYPE_INDEX]; // per allele non-log space probabilities of zero counts for this sample // for each allele calculate the total probability of genotypes containing at least one copy of the allele final double[] log10ProbabilityOfNonZeroAltAlleles = new double[numAlleles]; Arrays.fill(log10ProbabilityOfNonZeroAltAlleles, Double.NEGATIVE_INFINITY); for (int genotype = 0; genotype < glCalc.genotypeCount(); genotype++) { final double log10GenotypePosterior = log10GenotypePosteriors[genotype]; glCalc.genotypeAlleleCountsAt(genotype) .forEachAlleleIndexAndCount((alleleIndex, count) -> log10ProbabilityOfNonZeroAltAlleles[alleleIndex] = MathUtils .log10SumLog10(log10ProbabilityOfNonZeroAltAlleles[alleleIndex], log10GenotypePosterior)); } for (int allele = 0; allele < numAlleles; allele++) { // if prob of non hom ref == 1 up to numerical precision, short-circuit to avoid NaN if (log10ProbabilityOfNonZeroAltAlleles[allele] >= 0) { log10POfZeroCountsByAllele[allele] = Double.NEGATIVE_INFINITY; } else { log10POfZeroCountsByAllele[allele] += MathUtils .log10OneMinusPow10(log10ProbabilityOfNonZeroAltAlleles[allele]); } } } // unfortunately AFCalculationResult expects integers for the MLE. We really should emit the EM no-integer values // which are valuable (eg in CombineGVCFs) as the sufficient statistics of the Dirichlet posterior on allele frequencies final int[] integerAlleleCounts = Arrays.stream(alleleCounts).mapToInt(x -> (int) Math.round(x)).toArray(); final int[] integerAltAlleleCounts = Arrays.copyOfRange(integerAlleleCounts, 1, numAlleles); //skip the ref allele (index 0) final Map<Allele, Double> log10PRefByAllele = IntStream.range(1, numAlleles).boxed() .collect(Collectors.toMap(alleles::get, a -> log10POfZeroCountsByAllele[a])); // we compute posteriors here and don't have the same prior that AFCalculationResult expects. Therefore, we // give it our posterior as its "likelihood" along with a flat dummy prior final double[] dummyFlatPrior = { -1e-10, -1e-10 }; //TODO: HACK must be negative for AFCalcResult final double[] log10PosteriorOfNoVariantYesVariant = { log10PNoVariant, MathUtils.log10OneMinusPow10(log10PNoVariant) }; return new AFCalculationResult(integerAltAlleleCounts, DUMMY_N_EVALUATIONS, alleles, log10PosteriorOfNoVariantYesVariant, dummyFlatPrior, log10PRefByAllele); }
From source file:org.pmad.gmm.MyMixMNDEM.java
/** * Fit a mixture model to the data supplied to the constructor. * * The quality of the fit depends on the concavity of the data provided to * the constructor and the initial mixture provided to this function. If the * data has many local optima, multiple runs of the fitting function with * different initial mixtures may be required to find the optimal solution. * If a SingularMatrixException is encountered, it is possible that another * initialization would work.// w w w . jav a2 s .co m * * @param initialMixture Model containing initial values of weights and * multivariate normals * @param maxIterations Maximum iterations allowed for fit * @param threshold Convergence threshold computed as difference in * logLikelihoods between successive iterations * @throws SingularMatrixException if any component's covariance matrix is * singular during fitting * @throws NotStrictlyPositiveException if numComponents is less than one * or threshold is less than Double.MIN_VALUE * @throws DimensionMismatchException if initialMixture mean vector and data * number of columns are not equal */ public void fit(final MyMixMND initialMixture, final int maxIterations, final double threshold) throws SingularMatrixException, NotStrictlyPositiveException, DimensionMismatchException { if (maxIterations < 1) { throw new NotStrictlyPositiveException(maxIterations); } if (threshold < Double.MIN_VALUE) { throw new NotStrictlyPositiveException(threshold); } final int n = data.length; // Number of data columns. Jagged data already rejected in constructor, // so we can assume the lengths of each row are equal. final int numCols = data[0].length; final int k = initialMixture.getComponents().size(); final int numMeanColumns = initialMixture.getComponents().get(0).getSecond().getMeans().length; if (numMeanColumns != numCols) { throw new DimensionMismatchException(numMeanColumns, numCols); } int numIterations = 0; double previousLogLikelihood = 0d; logLikelihood = Double.NEGATIVE_INFINITY; // Initialize model to fit to initial mixture. fittedModel = new MyMixMND(initialMixture.getComponents()); while (numIterations++ <= maxIterations && FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { System.out.println(numIterations); previousLogLikelihood = logLikelihood; double sumLogLikelihood = 0d; // Mixture components final List<Pair<Double, MyMND>> components = fittedModel.getComponents(); // Weight and distribution of each component final double[] weights = new double[k]; final MyMND[] mvns = new MyMND[k]; for (int j = 0; j < k; j++) { weights[j] = components.get(j).getFirst(); mvns[j] = components.get(j).getSecond(); } // E-step: compute the data dependent parameters of the expectation // function. // The percentage of row's total density between a row and a // component final double[][] gamma = new double[n][k]; // Sum of gamma for each component final double[] gammaSums = new double[k]; // for (double d : weights) { // System.out.print(d+" "); // } // System.out.println(); // Sum of gamma times its row for each each component final double[][] gammaDataProdSums = new double[k][numCols]; for (int i = 0; i < n; i++) { final double rowDensity = fittedModel.density(data[i]); sumLogLikelihood += FastMath.log(rowDensity); for (int j = 0; j < k; j++) { gamma[i][j] = weights[j] * mvns[j].density(data[i]); if (rowDensity == 0) { if (gamma[i][j] == 0) { gamma[i][j] = weights[j]; } else { System.err.println("ele"); } } else { gamma[i][j] /= rowDensity; } gammaSums[j] += gamma[i][j]; for (int col = 0; col < numCols; col++) { gammaDataProdSums[j][col] += gamma[i][j] * data[i][col]; } } } if (Utils.hasNaN(gamma)) { System.out.println("gamma has NaN"); } logLikelihood = sumLogLikelihood / n; // M-step: compute the new parameters based on the expectation // function. final double[] newWeights = new double[k]; final double[][] newMeans = new double[k][numCols]; for (int j = 0; j < k; j++) { newWeights[j] = gammaSums[j] / n; for (int col = 0; col < numCols; col++) { newMeans[j][col] = gammaDataProdSums[j][col] / gammaSums[j]; } } // Compute new covariance matrices final RealMatrix[] newCovMats = new RealMatrix[k]; for (int j = 0; j < k; j++) { newCovMats[j] = new Array2DRowRealMatrix(numCols, numCols); } for (int i = 0; i < n; i++) { for (int j = 0; j < k; j++) { final RealMatrix vec = new Array2DRowRealMatrix(MathArrays.ebeSubtract(data[i], newMeans[j])); final RealMatrix dataCov = vec.multiply(vec.transpose()).scalarMultiply(gamma[i][j]); newCovMats[j] = newCovMats[j].add(dataCov); } } // Converting to arrays for use by fitted model final double[][][] newCovMatArrays = new double[k][numCols][numCols]; for (int j = 0; j < k; j++) { newCovMats[j] = newCovMats[j].scalarMultiply(1d / gammaSums[j]); newCovMatArrays[j] = newCovMats[j].getData(); // if (Utils.hasNaN(newCovMatArrays[j])) { // System.out.println("covmet nan"); // } } // if (Utils.hasNaN(newMeans)) { // System.out.println("neams nan"); // } // Update current model fittedModel = new MyMixMND(newWeights, newMeans, newCovMatArrays); } if (FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { // Did not converge before the maximum number of iterations // throw new ConvergenceException(); System.out.println("Did not converge"); } }