Java Code Examples for org.nd4j.linalg.api.ndarray.INDArray#div()

public static boolean checkDivManually(INDArray first, INDArray second, double maxRelativeDifference,
                double minAbsDifference) {
    //No apache commons element-wise division, but can do this manually

    INDArray result = first.div(second);
    long[] shape = first.shape();

    INDArray expected = Nd4j.zeros(first.shape());

    for (int i = 0; i < shape[0]; i++) {
        for (int j = 0; j < shape[1]; j++) {
            double v = first.getDouble(i, j) / second.getDouble(i, j);
            expected.putScalar(new int[] {i, j}, v);
        }
    }
    if (!checkShape(expected, result))
        return false;
    boolean ok = checkEntries(expected, result, maxRelativeDifference, minAbsDifference);
    if (!ok) {
        INDArray onCopies = Shape.toOffsetZeroCopy(first).mul(Shape.toOffsetZeroCopy(second));
        printFailureDetails(first, second, expected, result, onCopies, "div");
    }
    return ok;
}
 

@Test
public void testElementWiseOps() {
    INDArray n1 = Nd4j.scalar(1);
    INDArray n2 = Nd4j.scalar(2);
    INDArray nClone = n1.add(n2);
    assertEquals(Nd4j.scalar(3), nClone);
    INDArray n1PlusN2 = n1.add(n2);
    assertFalse(getFailureMessage(), n1PlusN2.equals(n1));

    INDArray n3 = Nd4j.scalar(3.0);
    INDArray n4 = Nd4j.scalar(4.0);
    INDArray subbed = n4.sub(n3);
    INDArray mulled = n4.mul(n3);
    INDArray div = n4.div(n3);

    assertFalse(subbed.equals(n4));
    assertFalse(mulled.equals(n4));
    assertEquals(Nd4j.scalar(1.0), subbed);
    assertEquals(Nd4j.scalar(12.0), mulled);
    assertEquals(Nd4j.scalar(1.333333333333333333333), div);
}
 

@Test
    public void testSimple(){
        Nd4j.create(1);
        for(DataType dt : new DataType[]{DataType.DOUBLE, DataType.FLOAT, DataType.HALF, DataType.INT, DataType.LONG}) {
//            System.out.println("----- " + dt + " -----");
            INDArray arr = Nd4j.ones(dt,1, 5);
//            System.out.println("Ones: " + arr);
            arr.assign(1.0);
//            System.out.println("assign(1.0): " + arr);
//            System.out.println("DIV: " + arr.div(8));
//            System.out.println("MUL: " + arr.mul(8));
//            System.out.println("SUB: " + arr.sub(8));
//            System.out.println("ADD: " + arr.add(8));
//            System.out.println("RDIV: " + arr.rdiv(8));
//            System.out.println("RSUB: " + arr.rsub(8));
            arr.div(8);
            arr.mul(8);
            arr.sub(8);
            arr.add(8);
            arr.rdiv(8);
            arr.rsub(8);
        }
    }
 

@Override
public INDArray computeGradient(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask) {
    if(!labels.equalShapes(preOutput)){
        Preconditions.throwEx("Labels and preOutput must have equal shapes: got shapes %s vs %s", labels.shape(), preOutput.shape());
    }
    labels = labels.castTo(preOutput.dataType());   //No-op if already correct dtype
    INDArray dLda = labels.div(labels.size(1));

    if (mask != null && LossUtil.isPerOutputMasking(dLda, mask)) {
        LossUtil.applyMask(labels, mask);
    }

    INDArray grad = activationFn.backprop(preOutput, dLda).getFirst();

    if (mask != null) {
        LossUtil.applyMask(grad, mask);
    }

    return grad;
}
 

public static boolean checkDivManually(INDArray first, INDArray second, double maxRelativeDifference,
                double minAbsDifference) {
    //No apache commons element-wise division, but can do this manually

    INDArray result = first.div(second);
    long[] shape = first.shape();

    INDArray expected = Nd4j.zeros(first.shape());

    for (int i = 0; i < shape[0]; i++) {
        for (int j = 0; j < shape[1]; j++) {
            double v = first.getDouble(i, j) / second.getDouble(i, j);
            expected.putScalar(new int[] {i, j}, v);
        }
    }
    if (!checkShape(expected, result))
        return false;
    boolean ok = checkEntries(expected, result, maxRelativeDifference, minAbsDifference);
    if (!ok) {
        INDArray onCopies = Shape.toOffsetZeroCopy(first).mul(Shape.toOffsetZeroCopy(second));
        printFailureDetails(first, second, expected, result, onCopies, "div");
    }
    return ok;
}
 

@Override
public INDArray computeGradient(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask) {
    if(!labels.equalShapes(preOutput)){
        Preconditions.throwEx("Labels and preOutput must have equal shapes: got shapes %s vs %s", labels.shape(), preOutput.shape());
    }
    labels = labels.castTo(preOutput.dataType());   //No-op if already correct dtype
    INDArray yHat = activationFn.getActivation(preOutput.dup(), true);
    INDArray yDivyhat = labels.div(yHat);
    INDArray dLda = yDivyhat.rsubi(1);

    if (mask != null && LossUtil.isPerOutputMasking(dLda, mask)) {
        //For *most* activation functions: we don't actually need to mask dL/da in addition to masking dL/dz later
        //but: some, like softmax, require both (due to dL/dz_i being a function of dL/da_j, for i != j)
        //We could add a special case for softmax (activationFn instanceof ActivationSoftmax) but that would be
        // error prone - though buy us a tiny bit of performance
        LossUtil.applyMask(dLda, mask);
    }

    INDArray gradients = activationFn.backprop(preOutput, dLda).getFirst(); //TODO activation functions with params

    if (mask != null) {
        LossUtil.applyMask(gradients, mask);
    }

    return gradients;
}
 

@Test
public void testElementWiseOps() {
    INDArray n1 = Nd4j.scalar(1);
    INDArray n2 = Nd4j.scalar(2);
    INDArray nClone = n1.add(n2);
    assertEquals(Nd4j.scalar(3), nClone);
    INDArray n1PlusN2 = n1.add(n2);
    assertFalse(getFailureMessage(), n1PlusN2.equals(n1));

    INDArray n3 = Nd4j.scalar(3);
    INDArray n4 = Nd4j.scalar(4);
    INDArray subbed = n4.sub(n3);
    INDArray mulled = n4.mul(n3);
    INDArray div = n4.div(n3);

    assertFalse(subbed.equals(n4));
    assertFalse(mulled.equals(n4));
    assertEquals(Nd4j.scalar(1), subbed);
    assertEquals(Nd4j.scalar(12), mulled);
    assertEquals(Nd4j.scalar(1.333333333333333333333), div);
}
 

@Test
public void testDivide() {
    INDArray two = Nd4j.create(new float[] {2, 2, 2, 2});
    INDArray div = two.div(two);
    assertEquals(getFailureMessage(), Nd4j.ones(DataType.FLOAT, 4), div);

    INDArray half = Nd4j.create(new float[] {0.5f, 0.5f, 0.5f, 0.5f}, new long[] {2, 2});
    INDArray divi = Nd4j.create(new float[] {0.3f, 0.6f, 0.9f, 0.1f}, new long[] {2, 2});
    INDArray assertion = Nd4j.create(new float[] {1.6666666f, 0.8333333f, 0.5555556f, 5}, new long[] {2, 2});
    INDArray result = half.div(divi);
    assertEquals(getFailureMessage(), assertion, result);
}
 

@Override
public INDArray computeGradient(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask) {
    if (labels.size(1) != preOutput.size(1)) {
        throw new IllegalArgumentException(
                        "Labels array numColumns (size(1) = " + labels.size(1) + ") does not match output layer"
                                        + " number of outputs (nOut = " + preOutput.size(1) + ") ");

    }
    INDArray yHat = activationFn.getActivation(preOutput.dup(), true);
    INDArray yDivyhat = labels.div(yHat);
    INDArray dLda = yDivyhat.rsubi(1);

    if (mask != null && LossUtil.isPerOutputMasking(dLda, mask)) {
        //For *most* activation functions: we don't actually need to mask dL/da in addition to masking dL/dz later
        //but: some, like softmax, require both (due to dL/dz_i being a function of dL/da_j, for i != j)
        //We could add a special case for softmax (activationFn instanceof ActivationSoftmax) but that would be
        // error prone - though buy us a tiny bit of performance
        LossUtil.applyMask(dLda, mask);
    }

    INDArray gradients = activationFn.backprop(preOutput, dLda).getFirst(); //TODO activation functions with params

    if (mask != null) {
        LossUtil.applyMask(gradients, mask);
    }

    return gradients;
}
 

/**
 * Returns the word vector divided by the norm2 of the array
 * @param word the word to get the matrix for
 * @return the looked up matrix
 */
public INDArray getWordVectorMatrixNormalized(String word) {
    INDArray r = getWordVectorMatrix(word);
    if (r == null)
        return null;

    return r.div(Nd4j.getBlasWrapper().nrm2(r));
}
 

@Test
public void testReverseFlow2() {
    CudaGridExecutioner executioner = ((CudaGridExecutioner) Nd4j.getExecutioner());

    INDArray n1 = Nd4j.scalar(1);
    INDArray n2 = Nd4j.scalar(2);
    INDArray n3 = Nd4j.scalar(3);
    INDArray n4 = Nd4j.scalar(4);

    System.out.println("0: ------------------------");

    INDArray nClone = n1.add(n2);

    assertEquals(Nd4j.scalar(3), nClone);
    INDArray n1PlusN2 = n1.add(n2);
    assertFalse(n1PlusN2.equals(n1));

    System.out.println("2: ------------------------");

    System.out.println(n4);

    INDArray subbed = n4.sub(n3);
    INDArray mulled = n4.mul(n3);
    INDArray div = n4.div(n3);

    System.out.println("Subbed: " + subbed);
    System.out.println("Mulled: " + mulled);
    System.out.println("Div: " + div);
    System.out.println("4: ------------------------");

    assertFalse(subbed.equals(n4));
    assertFalse(mulled.equals(n4));

    assertEquals(0, executioner.getQueueLength());

    assertEquals(Nd4j.scalar(1), subbed);
    assertEquals(Nd4j.scalar(12), mulled);
    assertEquals(Nd4j.scalar(1.333333333333333333333), div);
}
 

/**
 * Solves a linear system using Apache commons methods [mat].[x] = [vec]
 *
 * @param mat the coefficients matrix (must be square and full-rank)
 * @param vec the right hand side vector
 * @param singularityThreshold a threshold for detecting singularity
 * @return solution of the linear system
 */
public static INDArray linsolve(@Nonnull final INDArray mat, @Nonnull final INDArray vec,
                                final double singularityThreshold) {
    if (mat.isScalar()) {
        return vec.div(mat.getDouble(0));
    }
    if (!mat.isSquare()) {
        throw new IllegalArgumentException("invalid array: must be a square matrix");
    }
    final RealVector sol = new LUDecomposition(Nd4jApacheAdapterUtils.convertINDArrayToApacheMatrix(mat),
            singularityThreshold).getSolver().solve(Nd4jApacheAdapterUtils.convertINDArrayToApacheVector(vec));
    return Nd4j.create(sol.toArray(), vec.shape());
}
 

@Test
public void testPoolingEdgeCases(){
    //Average pooling with same mode: should we include the padded values, when deciding what to divide by?
    ///*** Note: Mode 2 is the "DL4J always divide by kH*kW" approach ***

    /*
    Input:
    [ 1, 2, 3
      4, 5, 6
      7, 8, 9 ]


     Kernel 2, stride 1
     outH = 3, outW = 3 (i.e., ceil(in/stride)
     totalHPad = (outH-1) * strideH + kH - inH = (3-1)*1 + 2 - 3 = 1
     topPad = 0, bottomPad = 1
     leftPad = 0, rightPad = 1
     */

    for( char inputOrder : new char[]{'c', 'f'}) {
        for( char outputOrder : new char[]{'c', 'f'}) {

            INDArray input = Nd4j.create(1, 1, 3, 3);
            input.get(point(0), point(0), all(), all())
                    .assign(Nd4j.linspace(1, 9, 9).reshape('c', 3, 3))
                    .dup(inputOrder);

            input = input.dup('c');

            INDArray input2 = Nd4j.create(new double[]{1,2,3,4,5,6,7,8,9}, new int[]{1,1,3,3}, 'c');//.dup(inputOrder);
            assertEquals(input, input2);

            input = input2;

            for( int i=0; i<3; i++){
                for( int j=0; j<3; j++ ){
                    System.out.print(input.getDouble(0,0,i,j) + ",");
                }
                System.out.println();
            }
            System.out.println();

            INDArray sums = Nd4j.create(new double[][]{
                    {(1 + 2 + 4 + 5), (2 + 3 + 5 + 6), (3 + 6)},
                    {(4 + 5 + 7 + 8), (5 + 6 + 8 + 9), (6 + 9)},
                    {(7 + 8), (8 + 9), (9)}
            });

            INDArray divEnabled = Nd4j.create(new double[][]{
                    {4, 4, 2},
                    {4, 4, 2},
                    {2, 2, 1}
            });

            INDArray expEnabled = sums.div(divEnabled);
            INDArray expDl4j = sums.div(4);

            //https://github.com/deeplearning4j/libnd4j/blob/master/include/ops/declarable/generic/convo/pooling/avgpool2d.cpp
            DynamicCustomOp op1 = DynamicCustomOp.builder("avgpool2d")
                    .addIntegerArguments(new int[]{2, 2, 1, 1, 0, 0, 1, 1, 1, 0, 0})   //ky, kx, sy, sx, py, px, dy, dx, isSameMode, ???, divisor, nchw
                    .addInputs(input)
                    .addOutputs(Nd4j.create(new int[]{1, 1, 3, 3}, outputOrder))
                    .build();

            DynamicCustomOp op2 = DynamicCustomOp.builder("avgpool2d")
                    .addIntegerArguments(new int[]{2, 2, 1, 1, 0, 0, 1, 1, 1, 2, 0})   //ky, kx, sy, sx, py, px, dy, dx, isSameMode, ???, divisor, nchw
                    .addInputs(input)
                    .addOutputs(Nd4j.create(new int[]{1, 1, 3, 3}, outputOrder))
                    .build();

            Nd4j.getExecutioner().exec(op1);
            Nd4j.getExecutioner().exec(op2);
            INDArray actEnabled = op1.getOutputArgument(0);
            INDArray actDl4j = op2.getOutputArgument(0);


            String msg = "inOrder=" + inputOrder + ", outOrder=" + outputOrder;
            assertEquals(msg, expDl4j, actDl4j.get(point(0), point(0), all(), all()));
            assertEquals(msg, expEnabled, actEnabled.get(point(0), point(0), all(), all()));
        }
    }
}
 

@Test
public void testForMultiDataSet() {
    DataSet dataSetA = knownDistVariedDataSet(new float[] {0.8f, 0.1f, 0.2f}, false);
    DataSet dataSetB = knownDistVariedDataSet(new float[] {0.2f, 0.9f, 0.8f}, true);

    HashMap<Integer, Double> targetDists = new HashMap<>();
    targetDists.put(0, 0.5); //balance inputA
    targetDists.put(1, 0.3); //inputB dist = 0.2%
    UnderSamplingByMaskingMultiDataSetPreProcessor maskingMultiDataSetPreProcessor =
                    new UnderSamplingByMaskingMultiDataSetPreProcessor(targetDists, window);
    maskingMultiDataSetPreProcessor.overrideMinorityDefault(1);

    MultiDataSet multiDataSet = fromDataSet(dataSetA, dataSetB);
    maskingMultiDataSetPreProcessor.preProcess(multiDataSet);

    INDArray labels;
    INDArray minorityCount;
    INDArray seqCount;
    INDArray minorityDist;
    //datasetA
    labels = multiDataSet.getLabels(0).reshape(minibatchSize * 2, longSeq).mul(multiDataSet.getLabelsMaskArray(0));
    minorityCount = labels.sum(1);
    seqCount = multiDataSet.getLabelsMaskArray(0).sum(1);
    minorityDist = minorityCount.div(seqCount);
    assertEquals(minorityDist.getDouble(1, 0), 0.5, tolerancePerc);
    assertEquals(minorityDist.getDouble(2, 0), 0.5, tolerancePerc);
    assertEquals(minorityDist.getDouble(4, 0), 0.5, tolerancePerc);
    assertEquals(minorityDist.getDouble(5, 0), 0.5, tolerancePerc);

    //datasetB - override is switched so grab index=0
    labels = multiDataSet.getLabels(1).get(NDArrayIndex.all(), NDArrayIndex.point(0), NDArrayIndex.all())
                    .mul(multiDataSet.getLabelsMaskArray(1));
    minorityCount = labels.sum(1);
    seqCount = multiDataSet.getLabelsMaskArray(1).sum(1);
    minorityDist = minorityCount.div(seqCount);
    assertEquals(minorityDist.getDouble(1, 0), 0.3, tolerancePerc);
    assertEquals(minorityDist.getDouble(2, 0), 0.3, tolerancePerc);
    assertEquals(minorityDist.getDouble(4, 0), 0.3, tolerancePerc);
    assertEquals(minorityDist.getDouble(5, 0), 0.3, tolerancePerc);

}
 

/**
 * This method returns the gradient of the cost function with respect to the
 * output from the previous layer.  For this cost function, the gradient
 * is derived from Bishop's paper "Mixture Density Networks" (1994) which
 * gives an elegant closed-form expression for the derivatives with respect
 * to each of the output components.
 * @param labels Labels to train on.
 * @param preOutput Output of neural network before applying the final activation function.
 * @param activationFn Activation function of output layer.
 * @param mask Mask to apply to gradients.
 * @return Gradient of cost function with respect to preOutput parameters.
 */
@Override
public INDArray computeGradient(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask) {
    long nSamples = labels.size(0);

    INDArray output = activationFn.getActivation(preOutput.dup(), false);

    MixtureDensityComponents mdc = extractComponents(output);

    INDArray gradient = Nd4j.zeros(nSamples, preOutput.columns());

    INDArray labelsMinusMu = labelsMinusMu(labels, mdc.mu);
    INDArray labelsMinusMuSquared = labelsMinusMu.mul(labelsMinusMu).sum(2);

    // This computes pi_i, see Bishop equation (30).
    // See http://www.plsyard.com/dealing-overflow-and-underflow-in-softmax-function/
    // this post for why we calculate the pi_i in this way.
    // With the exponential function here, we have to be very careful
    // about overflow/underflow considerations even with
    // fairly intermediate values.  Subtracting the max
    // here helps to ensure over/underflow does not happen here.
    // This isn't exactly a softmax because there's an 'alpha' coefficient
    // here, but the technique works, nonetheless.
    INDArray variance = mdc.sigma.mul(mdc.sigma);
    INDArray minustwovariance = variance.mul(2).negi();
    INDArray normalPart = mdc.alpha.div(Transforms.pow(mdc.sigma.mul(SQRT_TWO_PI), mLabelWidth));
    INDArray exponent = labelsMinusMuSquared.div(minustwovariance);
    INDArray exponentMax = exponent.max(1);
    exponent.subiColumnVector(exponentMax);
    INDArray pi = Transforms.exp(exponent).muli(normalPart);
    INDArray piDivisor = pi.sum(1);
    pi.diviColumnVector(piDivisor);

    // See Bishop equation (35)
    //INDArray dLdZAlpha = Nd4j.zeros(nSamples, nLabelsPerSample, mMixturesPerLabel); //mdc.alpha.sub(pi);
    INDArray dLdZAlpha = mdc.alpha.sub(pi);
    // See Bishop equation (38)
    INDArray dLdZSigma = (labelsMinusMuSquared.div(variance).subi(mLabelWidth)).muli(-1).muli(pi);
    // See Bishop equation (39)

    // This turned out to be way less efficient than
    // the simple 'for' loop here.
    //INDArray dLdZMu = pi
    //        .div(variance)
    //        .reshape(nSamples, mMixtures, 1)
    //        .repeat(2, mLabelWidth)
    //        .muli(labelsMinusMu)
    //        .negi()
    //        .reshape(nSamples, mMixtures * mLabelWidth);

    INDArray dLdZMu = Nd4j.create(nSamples, mMixtures, mLabelWidth);
    for (int k = 0; k < mLabelWidth; k++) {
        dLdZMu.put(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.all(), NDArrayIndex.point(k)},
                        labelsMinusMu.get(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.all(),
                                        NDArrayIndex.point(k)}).muli(pi).divi(variance).negi());
    }
    dLdZMu = dLdZMu.reshape(nSamples, mMixtures * mLabelWidth);

    // Place components of gradient into gradient holder.
    gradient.put(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.interval(0, mMixtures)}, dLdZAlpha);
    gradient.put(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.interval(mMixtures, mMixtures * 2)},
                    dLdZSigma);
    gradient.put(new INDArrayIndex[] {NDArrayIndex.all(),
                    NDArrayIndex.interval(mMixtures * 2, (mLabelWidth + 2) * mMixtures)}, dLdZMu);

    INDArray gradients = activationFn.backprop(preOutput, gradient).getFirst();

    if (mask != null) {
        LossUtil.applyMask(gradients, mask);
    }

    return gradients;
}
 

public INDArray calculate(INDArray X, int targetDimensions, double perplexity) {
    // pca hook
    if (usePca) {
        X = PCA.pca(X, Math.min(50, X.columns()), normalize);
    } else if (normalize) {
        X.subi(X.min(Integer.MAX_VALUE));
        X = X.divi(X.max(Integer.MAX_VALUE));
        X = X.subiRowVector(X.mean(0));
    }


    int n = X.rows();
    // FIXME: this is wrong, another distribution required here
    Y = Nd4j.randn(X.dataType(), X.rows(), targetDimensions);
    INDArray dY = Nd4j.zeros(n, targetDimensions);
    INDArray iY = Nd4j.zeros(n, targetDimensions);
    INDArray gains = Nd4j.ones(n, targetDimensions);

    boolean stopLying = false;
    logger.debug("Y:Shape is = " + Arrays.toString(Y.shape()));

    // compute P-values
    INDArray P = x2p(X, tolerance, perplexity);

    // do training
    for (int i = 0; i < maxIter; i++) {
        INDArray sumY = pow(Y, 2).sum(1).transpose();

        //Student-t distribution
        //also un normalized q
        // also known as num in original implementation
        INDArray qu = Y.mmul(Y.transpose()).muli(-2).addiRowVector(sumY).transpose().addiRowVector(sumY).addi(1)
                        .rdivi(1);

        //          doAlongDiagonal(qu,new Zero());

        INDArray Q = qu.div(qu.sumNumber().doubleValue());
        BooleanIndexing.replaceWhere(Q, 1e-12, Conditions.lessThan(1e-12));

        INDArray PQ = P.sub(Q).muli(qu);

        logger.debug("PQ shape is: " + Arrays.toString(PQ.shape()));
        logger.debug("PQ.sum(1) shape is: " + Arrays.toString(PQ.sum(1).shape()));

        dY = diag(PQ.sum(1)).subi(PQ).mmul(Y).muli(4);


        if (i < switchMomentumIteration) {
            momentum = initialMomentum;
        } else {
            momentum = finalMomentum;
        }

        gains = gains.add(.2).muli(dY.cond(Conditions.greaterThan(0)).neq(iY.cond(Conditions.greaterThan(0))))
                        .addi(gains.mul(0.8).muli(dY.cond(Conditions.greaterThan(0))
                                        .eq(iY.cond(Conditions.greaterThan(0)))));

        BooleanIndexing.replaceWhere(gains, minGain, Conditions.lessThan(minGain));

        INDArray gradChange = gains.mul(dY);

        gradChange.muli(learningRate);

        iY.muli(momentum).subi(gradChange);

        double cost = P.mul(log(P.div(Q), false)).sumNumber().doubleValue();
        logger.info("Iteration [" + i + "] error is: [" + cost + "]");

        Y.addi(iY);
        //  Y.addi(iY).subiRowVector(Y.mean(0));
        INDArray tiled = Nd4j.tile(Y.mean(0), new int[] {Y.rows(), 1});
        Y.subi(tiled);

        if (!stopLying && (i > maxIter / 2 || i >= stopLyingIteration)) {
            P.divi(4);
            stopLying = true;
        }
    }
    return Y;
}
 

@Test
public void testForMultiDataSet() {
    DataSet dataSetA = knownDistVariedDataSet(new float[] {0.8f, 0.1f, 0.2f}, false);
    DataSet dataSetB = knownDistVariedDataSet(new float[] {0.2f, 0.9f, 0.8f}, true);

    HashMap<Integer, Double> targetDists = new HashMap<>();
    targetDists.put(0, 0.5); //balance inputA
    targetDists.put(1, 0.3); //inputB dist = 0.2%
    UnderSamplingByMaskingMultiDataSetPreProcessor maskingMultiDataSetPreProcessor =
                    new UnderSamplingByMaskingMultiDataSetPreProcessor(targetDists, window);
    maskingMultiDataSetPreProcessor.overrideMinorityDefault(1);

    MultiDataSet multiDataSet = fromDataSet(dataSetA, dataSetB);
    maskingMultiDataSetPreProcessor.preProcess(multiDataSet);

    INDArray labels;
    INDArray minorityCount;
    INDArray seqCount;
    INDArray minorityDist;
    //datasetA
    labels = multiDataSet.getLabels(0).reshape(minibatchSize * 2, longSeq).mul(multiDataSet.getLabelsMaskArray(0));
    minorityCount = labels.sum(1);
    seqCount = multiDataSet.getLabelsMaskArray(0).sum(1);
    minorityDist = minorityCount.div(seqCount);
    assertEquals(minorityDist.getDouble(1), 0.5, tolerancePerc);
    assertEquals(minorityDist.getDouble(2), 0.5, tolerancePerc);
    assertEquals(minorityDist.getDouble(4), 0.5, tolerancePerc);
    assertEquals(minorityDist.getDouble(5), 0.5, tolerancePerc);

    //datasetB - override is switched so grab index=0
    labels = multiDataSet.getLabels(1).get(NDArrayIndex.all(), NDArrayIndex.point(0), NDArrayIndex.all())
                    .mul(multiDataSet.getLabelsMaskArray(1));
    minorityCount = labels.sum(1);
    seqCount = multiDataSet.getLabelsMaskArray(1).sum(1);
    minorityDist = minorityCount.div(seqCount);
    assertEquals(minorityDist.getDouble(1), 0.3, tolerancePerc);
    assertEquals(minorityDist.getDouble(2), 0.3, tolerancePerc);
    assertEquals(minorityDist.getDouble(4), 0.3, tolerancePerc);
    assertEquals(minorityDist.getDouble(5), 0.3, tolerancePerc);

}
 

/**
 * This method returns the gradient of the cost function with respect to the
 * output from the previous layer.  For this cost function, the gradient
 * is derived from Bishop's paper "Mixture Density Networks" (1994) which
 * gives an elegant closed-form expression for the derivatives with respect
 * to each of the output components.
 * @param labels Labels to train on.
 * @param preOutput Output of neural network before applying the final activation function.
 * @param activationFn Activation function of output layer.
 * @param mask Mask to apply to gradients.
 * @return Gradient of cost function with respect to preOutput parameters.
 */
@Override
public INDArray computeGradient(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask) {
    labels = labels.castTo(preOutput.dataType());   //No-op if already correct dtype
    long nSamples = labels.size(0);

    INDArray output = activationFn.getActivation(preOutput.dup(), false);

    MixtureDensityComponents mdc = extractComponents(output);

    INDArray gradient = Nd4j.zeros(nSamples, preOutput.columns());

    INDArray labelsMinusMu = labelsMinusMu(labels, mdc.mu);
    INDArray labelsMinusMuSquared = labelsMinusMu.mul(labelsMinusMu).sum(2);

    // This computes pi_i, see Bishop equation (30).
    // See http://www.plsyard.com/dealing-overflow-and-underflow-in-softmax-function/
    // this post for why we calculate the pi_i in this way.
    // With the exponential function here, we have to be very careful
    // about overflow/underflow considerations even with
    // fairly intermediate values.  Subtracting the max
    // here helps to ensure over/underflow does not happen here.
    // This isn't exactly a softmax because there's an 'alpha' coefficient
    // here, but the technique works, nonetheless.
    INDArray variance = mdc.sigma.mul(mdc.sigma);
    INDArray minustwovariance = variance.mul(2).negi();
    INDArray normalPart = mdc.alpha.div(Transforms.pow(mdc.sigma.mul(SQRT_TWO_PI), mLabelWidth));
    INDArray exponent = labelsMinusMuSquared.div(minustwovariance);
    INDArray exponentMax = exponent.max(1);
    exponent.subiColumnVector(exponentMax);
    INDArray pi = Transforms.exp(exponent).muli(normalPart);
    INDArray piDivisor = pi.sum(true,1);
    pi.diviColumnVector(piDivisor);

    // See Bishop equation (35)
    //INDArray dLdZAlpha = Nd4j.zeros(nSamples, nLabelsPerSample, mMixturesPerLabel); //mdc.alpha.sub(pi);
    INDArray dLdZAlpha = mdc.alpha.sub(pi);
    // See Bishop equation (38)
    INDArray dLdZSigma = (labelsMinusMuSquared.div(variance).subi(mLabelWidth)).muli(-1).muli(pi);
    // See Bishop equation (39)

    // This turned out to be way less efficient than
    // the simple 'for' loop here.
    //INDArray dLdZMu = pi
    //        .div(variance)
    //        .reshape(nSamples, mMixtures, 1)
    //        .repeat(2, mLabelWidth)
    //        .muli(labelsMinusMu)
    //        .negi()
    //        .reshape(nSamples, mMixtures * mLabelWidth);

    INDArray dLdZMu = Nd4j.create(nSamples, mMixtures, mLabelWidth);
    for (int k = 0; k < mLabelWidth; k++) {
        dLdZMu.put(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.all(), NDArrayIndex.point(k)},
                        labelsMinusMu.get(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.all(),
                                        NDArrayIndex.point(k)}).muli(pi).divi(variance).negi());
    }
    dLdZMu = dLdZMu.reshape(nSamples, mMixtures * mLabelWidth);

    // Place components of gradient into gradient holder.
    gradient.put(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.interval(0, mMixtures)}, dLdZAlpha);
    gradient.put(new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.interval(mMixtures, mMixtures * 2)},
                    dLdZSigma);
    gradient.put(new INDArrayIndex[] {NDArrayIndex.all(),
                    NDArrayIndex.interval(mMixtures * 2, (mLabelWidth + 2) * mMixtures)}, dLdZMu);

    INDArray gradients = activationFn.backprop(preOutput, gradient).getFirst();

    if (mask != null) {
        LossUtil.applyMask(gradients, mask);
    }

    return gradients;
}
 

/**
 * E-step update of read depth ($d_s$)
 *
 * @return a {@link SubroutineSignal} containing the update size (key: "error_norm")
 */
@EvaluatesRDD @UpdatesRDD @CachesRDD
public SubroutineSignal updateReadDepthPosteriorExpectations(final double admixingRatio,
                                                             final boolean neglectBiasCovariates) {
    mapWorkers(cb -> cb.cloneWithUpdatedCachesByTag(CoverageModelEMComputeBlock.CoverageModelICGCacheTag.E_STEP_D));
    cacheWorkers("after E-step for read depth initialization");
    /* map each compute block to their respective read depth estimation data (a triple of rank-1 INDArray's each
     * with STATE elements) and reduce by pairwise addition */
    final ImmutablePair<INDArray, INDArray> factors = mapWorkersAndReduce(
            cb -> cb.getReadDepthLatentPosteriorData(neglectBiasCovariates),
            (p1, p2) -> ImmutablePair.of(p1.left.add(p2.left), p1.right.add(p2.right)));

    /* put together */
    final INDArray numerator = factors.left;
    final INDArray denominator = factors.right;

    final INDArray newSampleMeanLogReadDepths = numerator.div(denominator);
    final INDArray newSampleVarLogReadDepths = Nd4j.ones(denominator.shape()).div(denominator);

    final INDArray newSampleMeanLogReadDepthsAdmixed;
    final INDArray newSampleVarLogReadDepthsAdmixed;

    /* admix */
    newSampleMeanLogReadDepthsAdmixed = newSampleMeanLogReadDepths
            .mul(admixingRatio)
            .addi(sampleMeanLogReadDepths.mul(1.0 - admixingRatio));
    newSampleVarLogReadDepthsAdmixed = newSampleVarLogReadDepths
            .mul(admixingRatio)
            .addi(sampleVarLogReadDepths.mul(1.0 - admixingRatio));

    /* calculate the error only using the change in the mean */
    final double errorNormInfinity = CoverageModelEMWorkspaceMathUtils.getINDArrayNormInfinity(
            newSampleMeanLogReadDepthsAdmixed.sub(sampleMeanLogReadDepths));

    /* update local copies of E[log(d_s)] and var[log(d_s)] */
    sampleMeanLogReadDepths.assign(newSampleMeanLogReadDepthsAdmixed);
    sampleVarLogReadDepths.assign(newSampleVarLogReadDepthsAdmixed);

    /* push E[log(d_s)] and var[log(d_s)] to all workers; they all need a copy */
    pushToWorkers(ImmutablePair.of(newSampleMeanLogReadDepths, newSampleVarLogReadDepths),
            (dat, cb) -> cb
                    .cloneWithUpdatedPrimitive(CoverageModelEMComputeBlock.CoverageModelICGCacheNode.log_d_s, dat.left)
                    .cloneWithUpdatedPrimitive(CoverageModelEMComputeBlock.CoverageModelICGCacheNode.var_log_d_s, dat.right));

    return SubroutineSignal.builder()
            .put(StandardSubroutineSignals.RESIDUAL_ERROR_NORM, errorNormInfinity)
            .build();
}
 

/**
 * Non Maximum Suppression - greedily selects the boxes with high confidence. Keep the boxes that have overlap area
 * below the threshold and discards the others.
 *
 * original code:
 *  - https://github.com/kpzhang93/MTCNN_face_detection_alignment/blob/master/code/codes/MTCNNv2/nms.m
 *  - https://github.com/davidsandberg/facenet/blob/master/src/align/detect_face.py#L687
 *
 * @param boxes nd array with bounding boxes: [[x1, y1, x2, y2 score]]
 * @param threshold NMS threshold -  retain overlap <= thresh
 * @param nmsType NMS method to apply. Available values ('Min', 'Union')
 * @return Returns the NMS result
 */
public static INDArray nonMaxSuppression(INDArray boxes, double threshold, NonMaxSuppressionType nmsType) {

	if (boxes.isEmpty()) {
		return Nd4j.empty();
	}

	// TODO Try to prevent following duplications!
	INDArray x1 = boxes.get(all(), point(0)).dup();
	INDArray y1 = boxes.get(all(), point(1)).dup();
	INDArray x2 = boxes.get(all(), point(2)).dup();
	INDArray y2 = boxes.get(all(), point(3)).dup();
	INDArray s = boxes.get(all(), point(4)).dup();

	//area = (x2 - x1 + 1) * (y2 - y1 + 1)
	INDArray area = (x2.sub(x1).add(1)).mul(y2.sub(y1).add(1));

	// sorted_s = np.argsort(s)
	INDArray sortedS = Nd4j.sortWithIndices(s, 0, SORT_ASCENDING)[0];

	INDArray pick = Nd4j.zerosLike(s);
	int counter = 0;

	while (sortedS.size(0) > 0) {

		if (sortedS.size(0) == 1) {
			pick.put(counter++, sortedS.dup());
			break;
		}

		long lastIndex = sortedS.size(0) - 1;
		INDArray i = sortedS.get(point(lastIndex), all()); // last element
		INDArray idx = sortedS.get(interval(0, lastIndex), all()).transpose(); // all until last excluding
		pick.put(counter++, i.dup());

		INDArray xx1 = Transforms.max(x1.get(idx), x1.get(i).getInt(0));
		INDArray yy1 = Transforms.max(y1.get(idx), y1.get(i).getInt(0));
		INDArray xx2 = Transforms.min(x2.get(idx), x2.get(i).getInt(0));
		INDArray yy2 = Transforms.min(y2.get(idx), y2.get(i).getInt(0));

		// w = np.maximum(0.0, xx2 - xx1 + 1)
		// h = np.maximum(0.0, yy2 - yy1 + 1)
		// inter = w * h
		INDArray w = Transforms.max(xx2.sub(xx1).add(1), 0.0f);
		INDArray h = Transforms.max(yy2.sub(yy1).add(1), 0.0f);
		INDArray inter = w.mul(h);

		// if method is 'Min':
		//   o = inter / np.minimum(area[i], area[idx])
		// else:
		//   o = inter / (area[i] + area[idx] - inter)
		int areaI = area.get(i).getInt(0);
		INDArray o = (nmsType == NonMaxSuppressionType.Min) ?
				inter.div(Transforms.min(area.get(idx), areaI)) :
				inter.div(area.get(idx).add(areaI).sub(inter));

		INDArray oIdx = MtcnnUtil.getIndexWhereVector(o, value -> value <= threshold);
		//INDArray oIdx = getIndexWhereVector2(o, Conditions.lessThanOrEqual(threshold));

		if (oIdx.isEmpty()) {
			break;
		}

		sortedS = Nd4j.expandDims(sortedS.get(oIdx), 0).transpose();
	}

	//pick = pick[0:counter]
	return (counter == 0) ? Nd4j.empty() : pick.get(interval(0, counter));
}