Java Code Examples for org.apache.commons.math3.distribution.UniformRealDistribution#sample()

public BM(int numClusters, int dimension, long randomSeed) {
    this.numClusters = numClusters;
    this.dimension = dimension;
    this.distributions = new BernoulliDistribution[numClusters][dimension];
    this.mixtureCoefficients = new double[numClusters];
    Arrays.fill(mixtureCoefficients,1.0/numClusters);
    this.logMixtureCoefficients = new double[numClusters];
    Arrays.fill(logMixtureCoefficients,Math.log(1.0/numClusters));
    Random random = new Random(randomSeed);
    RandomGenerator randomGenerator = RandomGeneratorFactory.createRandomGenerator(random);
    UniformRealDistribution uniform = new UniformRealDistribution(randomGenerator, 0.25,0.75);
    for (int k=0;k<numClusters;k++){
        for (int d=0;d<dimension;d++){
            double p = uniform.sample();
            distributions[k][d] = new BernoulliDistribution(p);
        }
    }
    this.logClusterConditioinalForEmpty = new double[numClusters];
    updateLogClusterConditioinalForEmpty();

    this.names = new ArrayList<>(dimension);
    for (int d=0;d<dimension;d++){
        names.add(""+d);
    }
}
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for ( int i = 0; i < sz; i++ ) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value( double x, double y ) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for ( int i = 0; i < xval.length; i++ ) {
            for ( int j = 0; j < yval.length; j++ ) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );

        final int numSamples = 50;
        final double tol = 5e-13;
        for ( int i = 0; i < numSamples; i++ ) {
            x = distX.sample();
            for ( int j = 0; j < numSamples; j++ ) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for ( int i = 0; i < sz; i++ ) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value( double x, double y ) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for ( int i = 0; i < xval.length; i++ ) {
            for ( int j = 0; j < yval.length; j++ ) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );

        final int numSamples = 50;
        final double tol = 5e-13;
        for ( int i = 0; i < numSamples; i++ ) {
            x = distX.sample();
            for ( int j = 0; j < numSamples; j++ ) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
 * @param numSamples Number of test samples.
 * @param tolerance Allowed tolerance on the interpolated value.
 * @param f Test function.
 * @param print Whether to print debugging output to the console.
 */
private void testInterpolation(int numSamples,
                               double tolerance,
                               BivariateFunction f,
                               boolean print) {
    final int sz = 21;
    final double[] xval = new double[sz];
    final double[] yval = new double[sz];
    // Coordinate values
    final double delta = 1d / (sz - 1);
    for (int i = 0; i < sz; i++) {
        xval[i] = -1 + 15 * i * delta;
        yval[i] = -20 + 30 * i * delta;
    }

    final double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    final BicubicInterpolator interpolator = new BicubicInterpolator();
    final BicubicInterpolatingFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;

    final RandomGenerator rng = new Well19937c();
    final UniformRealDistribution distX = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
    final UniformRealDistribution distY = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

    int count = 0;
    while (true) {
        x = distX.sample();
        y = distY.sample();
        if (!p.isValidPoint(x, y)) {
            if (print) {
                System.out.println("# " + x + " " + y);
            }
            continue;
        }

        if (count++ > numSamples) {
            break;
        }
        final double expected = f.value(x, y);
        final double actual = p.value(x, y);

        if (print) {
            System.out.println(x + " " + y + " " + expected + " " + actual);
        }

        Assert.assertEquals(expected, actual, tolerance);
    }
}
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }
        // Partial derivatives with respect to x
        double[][] dZdX = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdX[i][j] = 2;
            }
        }
        // Partial derivatives with respect to y
        double[][] dZdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdY[i][j] = -3;
            }
        }
        // Partial cross-derivatives
        double[][] dZdXdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdXdY[i][j] = 0;
            }
        }

        final BivariateFunction bcf
            = new BicubicSplineInterpolatingFunction(xval, yval, zval,
                                                     dZdX, dZdY, dZdXdY);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
                Assert.assertEquals(f.value(x, y),  bcf.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

AlleleFractionSimulatedData(final SampleLocatableMetadata metadata,
                            final AlleleFractionGlobalParameters globalParameters,
                            final int numSegments,
                            final double averageHetsPerSegment,
                            final double averageDepth,
                            final RandomGenerator rng) {
    final AlleleFractionState.MinorFractions minorFractions = new AlleleFractionState.MinorFractions(numSegments);
    final List<AllelicCount> allelicCounts = new ArrayList<>();
    final List<SimpleInterval> segments = new ArrayList<>();

    final PoissonDistribution segmentLengthGenerator = makePoisson(rng, averageHetsPerSegment);
    final PoissonDistribution readDepthGenerator = makePoisson(rng, averageDepth);
    final UniformRealDistribution minorFractionGenerator = new UniformRealDistribution(rng, 0.0, 0.5);

    final double meanBias = globalParameters.getMeanBias();
    final double biasVariance = globalParameters.getBiasVariance();
    final double outlierProbability = globalParameters.getOutlierProbability();

    //translate to ApacheCommons' parametrization of the gamma distribution
    final double gammaShape = meanBias * meanBias / biasVariance;
    final double gammaScale = biasVariance / meanBias;
    final GammaDistribution biasGenerator = new GammaDistribution(rng, gammaShape, gammaScale);

    //put each segment on its own chromosome and sort in sequence-dictionary order
    final List<String> chromosomes = IntStream.range(0, numSegments)
            .mapToObj(i -> metadata.getSequenceDictionary().getSequence(i).getSequenceName())
            .collect(Collectors.toList());

    for (final String chromosome : chromosomes) {
        // calculate the range of het indices for this segment
        final int numHetsInSegment = Math.max(MIN_HETS_PER_SEGMENT, segmentLengthGenerator.sample());

        final double minorFraction = minorFractionGenerator.sample();
        minorFractions.add(minorFraction);

        //we will put all the hets in this segment/chromosome at loci 1, 2, 3 etc
        segments.add(new SimpleInterval(chromosome, 1, numHetsInSegment));
        for (int het = 1; het < numHetsInSegment + 1; het++) {
            final double bias = biasGenerator.sample();

            //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
            final boolean isAltMinor = rng.nextDouble() < 0.5;
            final double altFraction =  isAltMinor ? minorFraction : 1 - minorFraction;

            //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
            final double pAlt;
            if (rng.nextDouble() < outlierProbability) {
                pAlt = rng.nextDouble();
            } else {
                pAlt = altFraction / (altFraction + (1 - altFraction) * bias);
            }

            final int numReads = readDepthGenerator.sample();
            final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
            final int numRefReads = numReads - numAltReads;
            allelicCounts.add(new AllelicCount(new SimpleInterval(chromosome, het, het), numRefReads, numAltReads));
        }
    }

    data = new AlleleFractionSegmentedData(
            new AllelicCountCollection(metadata, allelicCounts),
            new SimpleIntervalCollection(metadata, segments));
    trueState = new AlleleFractionState(meanBias, biasVariance, outlierProbability, minorFractions);
}
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }
        // Partial derivatives with respect to x
        double[][] dZdX = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdX[i][j] = 2;
            }
        }
        // Partial derivatives with respect to y
        double[][] dZdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdY[i][j] = -3;
            }
        }
        // Partial cross-derivatives
        double[][] dZdXdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdXdY[i][j] = 0;
            }
        }

        final BivariateFunction bcf
            = new BicubicSplineInterpolatingFunction(xval, yval, zval,
                                                     dZdX, dZdY, dZdXdY);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
                Assert.assertEquals(f.value(x, y),  bcf.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 251;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

public AlleleFractionSimulatedData(final double averageHetsPerSegment, final int numSegments,
        final double averageDepth, final double biasMean, final double biasVariance, final double outlierProbability) {
    rng.setSeed(RANDOM_SEED);
    this.numSegments = numSegments;
    final AlleleFractionState.MinorFractions minorFractions = new AlleleFractionState.MinorFractions(numSegments);
    final List<AllelicCount> alleleCounts = new ArrayList<>();
    final List<SimpleInterval> segments = new ArrayList<>();

    final PoissonDistribution segmentLengthGenerator = makePoisson(rng, averageHetsPerSegment);
    final PoissonDistribution readDepthGenerator = makePoisson(rng, averageDepth);
    final UniformRealDistribution minorFractionGenerator = new UniformRealDistribution(rng, 0.0, 0.5);

    //translate to ApacheCommons' parametrization of the gamma distribution
    final double gammaShape = biasMean * biasMean / biasVariance;
    final double gammaScale = biasVariance / biasMean;
    final GammaDistribution biasGenerator = new GammaDistribution(rng, gammaShape, gammaScale);

    //put each segment on its own chromosome and sort by lexicographical order
    final List<String> chromosomes = IntStream.range(0, numSegments).mapToObj(Integer::toString).collect(Collectors.toList());
    Collections.sort(chromosomes);

    for (final String chromosome : chromosomes) {
        // calculate the range of het indices for this segment
        final int numHetsInSegment = Math.max(MIN_HETS_PER_SEGMENT, segmentLengthGenerator.sample());

        final double minorFraction = minorFractionGenerator.sample();
        minorFractions.add(minorFraction);

        //we will put all the hets in this segment/chromosome at loci 1, 2, 3 etc
        segments.add(new SimpleInterval(chromosome, 1, numHetsInSegment + 1));
        for (int het = 1; het < numHetsInSegment + 1; het++) {
            final double bias = biasGenerator.sample();

            //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
            final boolean isAltMinor = rng.nextDouble() < 0.5;
            final double altFraction =  isAltMinor ? minorFraction : 1 - minorFraction;

            //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
            final double pAlt;
            if (rng.nextDouble() < outlierProbability) {
                truePhases.add(AlleleFractionIndicator.OUTLIER);
                pAlt = rng.nextDouble();
            } else {
                truePhases.add(isAltMinor ? AlleleFractionIndicator.ALT_MINOR : AlleleFractionIndicator.REF_MINOR);
                pAlt = altFraction / (altFraction + (1 - altFraction) * bias);
            }

            final int numReads = readDepthGenerator.sample();
            final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
            final int numRefReads = numReads - numAltReads;
            alleleCounts.add(new AllelicCount(new SimpleInterval(chromosome, het, het), numRefReads, numAltReads));
        }
    }

    final Genome genome = new Genome(TRIVIAL_TARGETS, alleleCounts);
    segmentedGenome = new SegmentedGenome(segments, genome);
    trueState = new AlleleFractionState(biasMean, biasVariance, outlierProbability, minorFractions);
}
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for ( int i = 0; i < sz; i++ ){
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value( double x, double y ) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for ( int i = 0; i < xval.length; i++ ) {
            for ( int j = 0; j < yval.length; j++ ) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );

        final int numSamples = 50;
        final double tol = 2e-14;
        for ( int i = 0; i < numSamples; i++ ) {
            x = distX.sample();
            for ( int j = 0; j < numSamples; j++ ) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }
        // Partial derivatives with respect to x
        double[][] dZdX = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdX[i][j] = 2;
            }
        }
        // Partial derivatives with respect to y
        double[][] dZdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdY[i][j] = -3;
            }
        }
        // Partial cross-derivatives
        double[][] dZdXdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdXdY[i][j] = 0;
            }
        }

        final BivariateFunction bcf
            = new BicubicSplineInterpolatingFunction(xval, yval, zval,
                                                     dZdX, dZdY, dZdXdY);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
                Assert.assertEquals(f.value(x, y),  bcf.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
 * @param minimumX Lower bound of interpolation range along the x-coordinate.
 * @param maximumX Higher bound of interpolation range along the x-coordinate.
 * @param minimumY Lower bound of interpolation range along the y-coordinate.
 * @param maximumY Higher bound of interpolation range along the y-coordinate.
 * @param numberOfElements Number of data points (along each dimension).
 * @param numberOfSamples Number of test points.
 * @param f Function to test.
 * @param meanTolerance Allowed average error (mean error on all interpolated values).
 * @param maxTolerance Allowed error on each interpolated value.
 */
private void testInterpolation(double minimumX,
                               double maximumX,
                               double minimumY,
                               double maximumY,
                               int numberOfElements,
                               int numberOfSamples,
                               BivariateFunction f,
                               double meanTolerance,
                               double maxTolerance) {
    double expected;
    double actual;
    double currentX;
    double currentY;
    final double deltaX = (maximumX - minimumX) / ((double) numberOfElements);
    final double deltaY = (maximumY - minimumY) / ((double) numberOfElements);
    final double[] xValues = new double[numberOfElements];
    final double[] yValues = new double[numberOfElements];
    final double[][] zValues = new double[numberOfElements][numberOfElements];

    for (int i = 0; i < numberOfElements; i++) {
        xValues[i] = minimumX + deltaX * (double) i;
        for (int j = 0; j < numberOfElements; j++) {
            yValues[j] = minimumY + deltaY * (double) j;
            zValues[i][j] = f.value(xValues[i], yValues[j]);
        }
    }

    final BivariateFunction interpolation
        = new PiecewiseBicubicSplineInterpolatingFunction(xValues,
                                                          yValues,
                                                          zValues);

    for (int i = 0; i < numberOfElements; i++) {
        currentX = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            currentY = yValues[j];
            expected = f.value(currentX, currentY);
            actual = interpolation.value(currentX, currentY);
            Assert.assertTrue(Precision.equals(expected, actual));
        }
    }

    final RandomGenerator rng = new Well19937c(1234567L);
    final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]);
    final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]);

    double sumError = 0;
    for (int i = 0; i < numberOfSamples; i++) {
        currentX = distX.sample();
        currentY = distY.sample();
        expected = f.value(currentX, currentY);
        actual = interpolation.value(currentX, currentY);
        sumError += FastMath.abs(actual - expected);
        Assert.assertEquals(expected, actual, maxTolerance);
    }

    final double meanError = sumError / numberOfSamples;
    Assert.assertEquals(0, meanError, meanTolerance);
}
 

CopyRatioSimulatedData(final SampleLocatableMetadata metadata,
                       final double variance,
                       final double outlierProbability,
                       final int numSegments,
                       final double averageIntervalsPerSegment,
                       final RandomGenerator rng) {
    final List<Double> segmentMeans = new ArrayList<>(numSegments);
    final List<Boolean> outlierIndicators = new ArrayList<>();
    final List<CopyRatio> copyRatios = new ArrayList<>();
    final List<SimpleInterval> segments = new ArrayList<>();

    final double standardDeviation = Math.sqrt(variance);

    final PoissonDistribution segmentLengthGenerator = makePoisson(rng, averageIntervalsPerSegment);
    final UniformRealDistribution segmentMeanGenerator = new UniformRealDistribution(rng, SEGMENT_MEAN_MIN, SEGMENT_MEAN_MAX);
    final UniformRealDistribution outlierGenerator = new UniformRealDistribution(rng, LOG2_COPY_RATIO_MIN, LOG2_COPY_RATIO_MAX);

    //put each segment on its own chromosome and sort in sequence-dictionary order
    final List<String> chromosomes = IntStream.range(0, numSegments)
            .mapToObj(i -> metadata.getSequenceDictionary().getSequence(i).getSequenceName())
            .collect(Collectors.toList());

    for (final String chromosome : chromosomes) {
        // calculate the range of interval indices for this segment
        final int numIntervalsInSegment = Math.max(MIN_INTERVALS_PER_SEGMENT, segmentLengthGenerator.sample());

        final double segmentMean = segmentMeanGenerator.sample();
        segmentMeans.add(segmentMean);

        //we will put all the intervals in this segment/chromosome at loci 1, 2, 3 etc
        segments.add(new SimpleInterval(chromosome, 1, numIntervalsInSegment));
        for (int interval = 1; interval < numIntervalsInSegment + 1; interval++) {

            final double log2CopyRatio;
            if (rng.nextDouble() < outlierProbability) {
                outlierIndicators.add(true);
                log2CopyRatio = outlierGenerator.sample();
            } else {
                outlierIndicators.add(false);
                log2CopyRatio = segmentMean + rng.nextGaussian() * standardDeviation;
            }

            copyRatios.add(new CopyRatio(new SimpleInterval(chromosome, interval, interval), log2CopyRatio));
        }
    }

    data = new CopyRatioSegmentedData(
            new CopyRatioCollection(metadata, copyRatios),
            new SimpleIntervalCollection(metadata, segments));
    trueState = new CopyRatioState(variance, outlierProbability, new CopyRatioState.SegmentMeans(segmentMeans), new CopyRatioState.OutlierIndicators(outlierIndicators));
}
 

/**
 * @param minimumX Lower bound of interpolation range along the x-coordinate.
 * @param maximumX Higher bound of interpolation range along the x-coordinate.
 * @param minimumY Lower bound of interpolation range along the y-coordinate.
 * @param maximumY Higher bound of interpolation range along the y-coordinate.
 * @param numberOfElements Number of data points (along each dimension).
 * @param numberOfSamples Number of test points.
 * @param f Function to test.
 * @param dfdx Partial derivative w.r.t. x of the function to test.
 * @param dfdy Partial derivative w.r.t. y of the function to test.
 * @param d2fdxdy Second partial cross-derivative of the function to test.
 * @param meanTolerance Allowed average error (mean error on all interpolated values).
 * @param maxTolerance Allowed error on each interpolated value.
 */
private void testInterpolation(double minimumX,
                               double maximumX,
                               double minimumY,
                               double maximumY,
                               int numberOfElements,
                               int numberOfSamples,
                               BivariateFunction f,
                               BivariateFunction dfdx,
                               BivariateFunction dfdy,
                               BivariateFunction d2fdxdy,
                               double meanTolerance,
                               double maxTolerance,
                               boolean print) {
    double expected;
    double actual;
    double currentX;
    double currentY;
    final double deltaX = (maximumX - minimumX) / numberOfElements;
    final double deltaY = (maximumY - minimumY) / numberOfElements;
    final double[] xValues = new double[numberOfElements];
    final double[] yValues = new double[numberOfElements];
    final double[][] zValues = new double[numberOfElements][numberOfElements];
    final double[][] dzdx = new double[numberOfElements][numberOfElements];
    final double[][] dzdy = new double[numberOfElements][numberOfElements];
    final double[][] d2zdxdy = new double[numberOfElements][numberOfElements];

    for (int i = 0; i < numberOfElements; i++) {
        xValues[i] = minimumX + deltaX * i;
        final double x = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            yValues[j] = minimumY + deltaY * j;
            final double y = yValues[j];
            zValues[i][j] = f.value(x, y);
            dzdx[i][j] = dfdx.value(x, y);
            dzdy[i][j] = dfdy.value(x, y);
            d2zdxdy[i][j] = d2fdxdy.value(x, y);
        }
    }

    final BivariateFunction interpolation
        = new BicubicInterpolatingFunction(xValues,
                                           yValues,
                                           zValues,
                                           dzdx,
                                           dzdy,
                                           d2zdxdy);

    for (int i = 0; i < numberOfElements; i++) {
        currentX = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            currentY = yValues[j];
            expected = f.value(currentX, currentY);
            actual = interpolation.value(currentX, currentY);
            Assert.assertTrue("On data point: " + expected + " != " + actual,
                              Precision.equals(expected, actual));
        }
    }

    final RandomGenerator rng = new Well19937c(1234567L);
    final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]);
    final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]);

    double sumError = 0;
    for (int i = 0; i < numberOfSamples; i++) {
        currentX = distX.sample();
        currentY = distY.sample();
        expected = f.value(currentX, currentY);

        if (print) {
            System.out.println(currentX + " " + currentY + " -> ");
        }

        actual = interpolation.value(currentX, currentY);
        sumError += FastMath.abs(actual - expected);

        if (print) {
            System.out.println(actual + " (diff=" + (expected - actual) + ")");
        }

        Assert.assertEquals(expected, actual, maxTolerance);
    }

    final double meanError = sumError / numberOfSamples;
    Assert.assertEquals(0, meanError, meanTolerance);
}
 

private void testInterpolation( double minimumX, double maximumX, int numberOfElements, int numberOfSamples,
                                UnivariateFunction f, double tolerance, double maxTolerance )
{
    double expected;
    double actual;
    double currentX;
    final double delta = ( maximumX - minimumX ) / ( (double) numberOfElements );
    double xValues[] = new double[numberOfElements];
    double yValues[] = new double[numberOfElements];

    for ( int i = 0; i < numberOfElements; i++ )
    {
        xValues[i] = minimumX + delta * (double) i;
        yValues[i] = f.value( xValues[i] );
    }

    UnivariateFunction interpolation = new AkimaSplineInterpolator().interpolate( xValues, yValues );

    for ( int i = 0; i < numberOfElements; i++ )
    {
        currentX = xValues[i];
        expected = f.value( currentX );
        actual = interpolation.value( currentX );
        assertTrue( Precision.equals( expected, actual ) );
    }

    final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed.
    final UniformRealDistribution distX =
        new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] );

    double sumError = 0;
    for ( int i = 0; i < numberOfSamples; i++ )
    {
        currentX = distX.sample();
        expected = f.value( currentX );
        actual = interpolation.value( currentX );
        sumError += FastMath.abs( actual - expected );
        assertEquals( expected, actual, maxTolerance );
    }

    assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance );
}
 

/**
 * @param minimumX Lower bound of interpolation range along the x-coordinate.
 * @param maximumX Higher bound of interpolation range along the x-coordinate.
 * @param minimumY Lower bound of interpolation range along the y-coordinate.
 * @param maximumY Higher bound of interpolation range along the y-coordinate.
 * @param numberOfElements Number of data points (along each dimension).
 * @param numberOfSamples Number of test points.
 * @param f Function to test.
 * @param dfdx Partial derivative w.r.t. x of the function to test.
 * @param dfdy Partial derivative w.r.t. y of the function to test.
 * @param d2fdxdy Second partial cross-derivative of the function to test.
 * @param meanTolerance Allowed average error (mean error on all interpolated values).
 * @param maxTolerance Allowed error on each interpolated value.
 */
private void testInterpolation(double minimumX,
                               double maximumX,
                               double minimumY,
                               double maximumY,
                               int numberOfElements,
                               int numberOfSamples,
                               BivariateFunction f,
                               BivariateFunction dfdx,
                               BivariateFunction dfdy,
                               BivariateFunction d2fdxdy,
                               double meanTolerance,
                               double maxTolerance,
                               boolean print) {
    double expected;
    double actual;
    double currentX;
    double currentY;
    final double deltaX = (maximumX - minimumX) / numberOfElements;
    final double deltaY = (maximumY - minimumY) / numberOfElements;
    final double[] xValues = new double[numberOfElements];
    final double[] yValues = new double[numberOfElements];
    final double[][] zValues = new double[numberOfElements][numberOfElements];
    final double[][] dzdx = new double[numberOfElements][numberOfElements];
    final double[][] dzdy = new double[numberOfElements][numberOfElements];
    final double[][] d2zdxdy = new double[numberOfElements][numberOfElements];

    for (int i = 0; i < numberOfElements; i++) {
        xValues[i] = minimumX + deltaX * i;
        final double x = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            yValues[j] = minimumY + deltaY * j;
            final double y = yValues[j];
            zValues[i][j] = f.value(x, y);
            dzdx[i][j] = dfdx.value(x, y);
            dzdy[i][j] = dfdy.value(x, y);
            d2zdxdy[i][j] = d2fdxdy.value(x, y);
        }
    }

    final BivariateFunction interpolation
        = new BicubicInterpolatingFunction(xValues,
                                           yValues,
                                           zValues,
                                           dzdx,
                                           dzdy,
                                           d2zdxdy);

    for (int i = 0; i < numberOfElements; i++) {
        currentX = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            currentY = yValues[j];
            expected = f.value(currentX, currentY);
            actual = interpolation.value(currentX, currentY);
            Assert.assertTrue("On data point: " + expected + " != " + actual,
                              Precision.equals(expected, actual));
        }
    }

    final RandomGenerator rng = new Well19937c(1234567L);
    final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]);
    final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]);

    double sumError = 0;
    for (int i = 0; i < numberOfSamples; i++) {
        currentX = distX.sample();
        currentY = distY.sample();
        expected = f.value(currentX, currentY);

        if (print) {
            System.out.println(currentX + " " + currentY + " -> ");
        }

        actual = interpolation.value(currentX, currentY);
        sumError += FastMath.abs(actual - expected);

        if (print) {
            System.out.println(actual + " (diff=" + (expected - actual) + ")");
        }

        Assert.assertEquals(expected, actual, maxTolerance);
    }

    final double meanError = sumError / numberOfSamples;
    Assert.assertEquals(0, meanError, meanTolerance);
}
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for ( int i = 0; i < sz; i++ ){
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value( double x, double y ) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for ( int i = 0; i < xval.length; i++ ) {
            for ( int j = 0; j < yval.length; j++ ) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );

        final int numSamples = 50;
        final double tol = 2e-14;
        for ( int i = 0; i < numSamples; i++ ) {
            x = distX.sample();
            for ( int j = 0; j < numSamples; j++ ) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }