Java Code Examples for org.nd4j.linalg.ops.transforms.Transforms#exp()

/**
 * Computes the gradient of the NLL for PL networks w.r.t. the k-th dyad according to equation (28) in [1].
 * @param plNetOutputs  The outputs for M_n dyads generated by a PLNet's output layer in order of their ranking (from best to worst).
 * @param k				The ranking position with respect to which the gradient should be computed. Assumes zero-based indices, unlike the paper.
 * @return				The gradient of the NLL loss w.r.t. the k-th dyad in the ranking.
 */
public static INDArray computeLossGradient(INDArray plNetOutputs, int k) {
	if (!(plNetOutputs.isRowVector()) || plNetOutputs.size(1) < 2 || k < 0 || k >= plNetOutputs.size(1)) {
		throw new IllegalArgumentException("Input has to be a row vector of 2 or more elements. And k has to be a valid index of plNetOutputs.");
	}
	long dyadRankingLength = plNetOutputs.size(1);
	double errorGradient = 0;
	for (int m = 0; m <= k; m++) {
		INDArray innerSumSlice = plNetOutputs.get(NDArrayIndex.interval(m, dyadRankingLength));
		innerSumSlice = Transforms.exp(innerSumSlice);
		double innerSum = innerSumSlice.sum(1).getDouble(0);
		errorGradient += Math.exp(plNetOutputs.getDouble(k)) / innerSum;
	}
	errorGradient -= 1;
	return Nd4j.create(new double[] {errorGradient});
}
 

@Test
public void testLogSumExp2() {

    for (int dim = 0; dim <= 2; dim++) {
        Nd4j.getRandom().setSeed(12345);
        INDArray inputArr = Nd4j.rand(DataType.DOUBLE, 3, 4, 5);
        SameDiff sd = SameDiff.create();
        SDVariable in = sd.var(inputArr);
        SDVariable lse = sd.math().logSumExp(in, dim);

        INDArray exp = Transforms.exp(inputArr, true);
        INDArray sum = exp.sum(dim);
        INDArray log = Transforms.log(sum);

        OpValidation.validate(new TestCase(sd)
                .expectedOutput(lse.name(), log)
                .gradientCheck(true));
    }
}
 

@Override
public INDArray op(INDArray neighbourhoodFeatures, INDArray nodeFeature) {
    double sigma = 16;
    final INDArray norm2 = Transforms.pow(neighbourhoodFeatures.subRowVector(nodeFeature), 2).sum(0);
    norm2.divi(-sigma * sigma);
    return Transforms.exp(norm2);
}
 

@Override
public INDArray generateRandom(INDArray preOutDistributionParams) {
    INDArray output = preOutDistributionParams.dup();
    activationFn.getActivation(output, true);

    val size = output.size(1) / 2;
    INDArray mean = output.get(NDArrayIndex.all(), NDArrayIndex.interval(0, size));
    INDArray logStdevSquared = output.get(NDArrayIndex.all(), NDArrayIndex.interval(size, 2 * size));

    INDArray sigma = Transforms.exp(logStdevSquared, true);
    Transforms.sqrt(sigma, false);

    INDArray e = Nd4j.randn(sigma.shape());
    return e.muli(sigma).addi(mean); //mu + sigma * N(0,1) ~ N(mu,sigma^2)
}
 

private INDArray[] calcLogProbArrayExConstants(INDArray x, INDArray preOutDistributionParams) {
    INDArray output = preOutDistributionParams.dup();
    activationFn.getActivation(output, false);

    val size = output.size(1) / 2;
    INDArray mean = output.get(NDArrayIndex.all(), NDArrayIndex.interval(0, size));
    INDArray logStdevSquared = output.get(NDArrayIndex.all(), NDArrayIndex.interval(size, 2 * size));

    INDArray sigmaSquared = Transforms.exp(logStdevSquared, true);
    INDArray lastTerm = x.sub(mean.castTo(x.dataType()));
    lastTerm.muli(lastTerm);
    lastTerm.divi(sigmaSquared.castTo(lastTerm.dataType())).divi(2);

    return new INDArray[] {logStdevSquared, lastTerm};
}
 

@Override
public INDArray generateAtMean(INDArray preOutDistributionParams) {
    //Input: gamma = log(lambda)    ->  lambda = exp(gamma)
    //Mean for exponential distribution: 1/lambda

    INDArray gamma = activationFn.getActivation(preOutDistributionParams.dup(), false);

    INDArray lambda = Transforms.exp(gamma, true);
    return lambda.rdivi(1.0); //mean = 1.0 / lambda
}
 

@Test
public void testLogSumExp() {
    Nd4j.getRandom().setSeed(12345);
    INDArray inputArr = Nd4j.rand(DataType.FLOAT, 1, 4);
    SameDiff sd = SameDiff.create();
    SDVariable in = sd.var(inputArr);
    SDVariable lse = sd.math().logSumExp(in);
    INDArray out = lse.eval();

    INDArray exp = Transforms.exp(inputArr, true);
    INDArray sum = exp.sum();
    INDArray log = Transforms.log(sum);
    assertEquals(log, out);
}
 

@Test
public void testExp() {
    INDArray n = Nd4j.create(new double[] {1, 2, 3, 4});
    INDArray assertion = Nd4j.create(new double[] {2.71828183f, 7.3890561f, 20.08553692f, 54.59815003f});
    INDArray exped = Transforms.exp(n);
    assertEquals(assertion, exped);
}
 

@Override
public INDArray gradient(INDArray x, INDArray preOutDistributionParams) {
    //p(x) = lambda * exp( -lambda * x)
    //logp(x) = log(lambda) - lambda * x = gamma - lambda * x
    //dlogp(x)/dgamma = 1 - lambda * x      (or negative of this for d(-logp(x))/dgamma

    INDArray gamma = activationFn.getActivation(preOutDistributionParams.dup(), true);

    INDArray lambda = Transforms.exp(gamma, true);
    INDArray dLdx = x.mul(lambda).subi(1.0);

    //dL/dz
    return activationFn.backprop(preOutDistributionParams.dup(), dLdx).getFirst();
}
 

private List<List<HiddenStateSegmentRecord<STATE, Target>>> getCopyRatioSegmentsLocal() {
    final List<List<CoverageModelCopyRatioEmissionData>> copyRatioEmissionData = fetchCopyRatioEmissionDataLocal();
    final INDArray sampleReadDepths = Transforms.exp(sampleMeanLogReadDepths, true);
    return sampleIndexStream()
            .mapToObj(si -> {
                final CopyRatioCallingMetadata metadata = CopyRatioCallingMetadata.builder()
                        .sampleName(processedSampleNameList.get(si))
                        .sampleSexGenotypeData(processedSampleSexGenotypeData.get(si))
                        .sampleCoverageDepth(sampleReadDepths.getDouble(si))
                        .emissionCalculationStrategy(EmissionCalculationStrategy.HYBRID_POISSON_GAUSSIAN)
                        .build();
                return copyRatioExpectationsCalculator.getCopyRatioHMMResults(metadata,
                        processedTargetList, copyRatioEmissionData.get(si));
            })
            /* segment each sample individually */
            .map(result -> {
                final HMMSegmentProcessor<CoverageModelCopyRatioEmissionData, STATE, Target> processor =
                        new HMMSegmentProcessor<>(
                                Collections.singletonList(result.getMetaData().getSampleName()),
                                Collections.singletonList(result.getMetaData().getSampleSexGenotypeData()),
                                referenceStateFactory,
                                Collections.singletonList(new HashedListTargetCollection<>(processedTargetList)),
                                Collections.singletonList(result.getForwardBackwardResult()),
                                Collections.singletonList(result.getViterbiResult()));
                return processor.getSegmentsAsList();
            })
            .collect(Collectors.toList());
}
 

/**
 * Computes the NLL for PL networks according to equation (27) in [1].
 * 
 * @param 	plNetOutputs The outputs for M_n dyads generated by a PLNet's output layer in order of their ranking (from best to worst).
 * @return	The NLL loss for the given PLNet outputs.
 */
public static INDArray computeLoss(INDArray plNetOutputs) {
	if (!(plNetOutputs.isRowVector()) || plNetOutputs.size(1) < 2 ) {
		throw new IllegalArgumentException("Input has to be a row vector of 2 or more elements.");
	}
	long dyadRankingLength = plNetOutputs.size(1);
	double loss = 0;
	for (int m = 0; m <= dyadRankingLength - 2; m++) {
		INDArray innerSumSlice = plNetOutputs.get(NDArrayIndex.interval(m, dyadRankingLength));
		innerSumSlice = Transforms.exp(innerSumSlice);
		loss += Transforms.log(innerSumSlice.sum(1)).getDouble(0);
	}
	loss -= plNetOutputs.get(NDArrayIndex.interval(0, dyadRankingLength - 1)).sum(1).getDouble(0);
	return Nd4j.create(new double[]{loss});
}
 

@Override
public INDArray ndOp(INDArray features, INDArray adjacencyMatrix) {
    double sigma = 16;
    INDArray[] sumsOfSquareDiffs = new INDArray[adjacencyMatrix.rows()];
    for (int node = 0; node < adjacencyMatrix.rows(); node++) {
        INDArray column = adjacencyMatrix.getColumn(node);
        INDArray repeat = features.getRow(node).repeat(0, features.rows()).muliColumnVector(column);
        INDArray sub = repeat.sub(features.mulColumnVector(column));
        sumsOfSquareDiffs[node] = Transforms.pow(sub, 2).sum(0);
    }
    INDArray sumOfSquareDiffs = Nd4j.vstack(sumsOfSquareDiffs).muli(-(1d / Math.pow(sigma, 2)));
    return Transforms.exp(sumOfSquareDiffs);
}
 

private void calculate(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask, INDArray scoreOutput, INDArray gradientOutput) {
    if (scoreOutput == null && gradientOutput == null) {
        throw new IllegalArgumentException("You have to provide at least one of scoreOutput or gradientOutput!");
    }
    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) + ") ");

    }
    final INDArray postOutput = activationFn.getActivation(preOutput.dup(), true);

    final INDArray positive = labels;
    final INDArray negative = labels.eq(0.0);
    final INDArray normFactor = negative.sum(1).muli(positive.sum(1));


    long examples = positive.size(0);
    for (int i = 0; i < examples; i++) {
        final INDArray locCfn = postOutput.getRow(i);
        final long[] shape = locCfn.shape();

        final INDArray locPositive = positive.getRow(i);
        final INDArray locNegative = negative.getRow(i);
        final Double locNormFactor = normFactor.getDouble(i);

        final INDArray operandA = Nd4j.ones(shape[1], shape[0]).mmul(locCfn);
        final INDArray operandB = operandA.transpose();

        final INDArray pairwiseSub = Transforms.exp(operandA.sub(operandB));

        final INDArray selection = locPositive.transpose().mmul(locNegative);

        final INDArray classificationDifferences = pairwiseSub.muli(selection).divi(locNormFactor);

        if (scoreOutput != null) {
            if (mask != null) {
                final INDArray perLabel = classificationDifferences.sum(0);
                LossUtil.applyMask(perLabel, mask.getRow(i));
                perLabel.sum(scoreOutput.getRow(i), 0);
            } else {
                classificationDifferences.sum(scoreOutput.getRow(i), 0, 1);
            }
        }

        if (gradientOutput != null) {
            gradientOutput.getRow(i).assign(classificationDifferences.sum(0).addi(classificationDifferences.sum(1).transposei().negi()));
        }
    }

    if (gradientOutput != null) {
        gradientOutput.assign(activationFn.backprop(preOutput.dup(), gradientOutput).getFirst());
        //multiply with masks, always
        if (mask != null) {
            LossUtil.applyMask(gradientOutput, mask);
        }
    }
}
 

public MixtureDensityComponents extractComponents(INDArray output) {
    long outputSize = output.size(1);
    if (outputSize != (mLabelWidth + 2) * mMixtures) {
        throw new IllegalArgumentException(
                        "Network output size " + outputSize + " must be (labels+2)*mixtures where labels = "
                                        + mLabelWidth + " and mixtures = " + mMixtures);
    }

    MixtureDensityComponents mdc = new MixtureDensityComponents();

    // Output is 2 dimensional (samples, labels)
    //
    // For each label vector of length 'labels', we will have
    // an output vector of length '(labels + 2) * nMixtures.
    // The first nMixtures outputs will correspond to the 'alpha' for each mixture.
    // The second nMixtures outputs will correspond to the 'sigma' and the last nMixtures*labels
    // will correspond to the 'mu' (mean) of the output.

    // Reorganize these.
    // alpha = samples, 0 to nMixtures
    // mu = samples, nMixtures to 2*nMixtures
    // sigma = samples, 2*nMixtures to (labels + 2)*nMixtures
    // Alpha is then sub-divided through reshape by mixtures per label and samples.

    mdc.alpha = output.get(NDArrayIndex.all(), NDArrayIndex.interval(0, mMixtures));
    mdc.sigma = output.get(NDArrayIndex.all(), NDArrayIndex.interval(mMixtures, 2 * mMixtures));
    mdc.mu = output.get(NDArrayIndex.all(), NDArrayIndex.interval(2 * mMixtures, (mLabelWidth + 2) * mMixtures))
                    .reshape(output.size(0), mMixtures, mLabelWidth);

    // Alpha is a softmax because
    // the alpha should all sum to 1 for a given gaussian mixture.
    mdc.alpha = Nd4j.getExecutioner().execAndReturn(new OldSoftMax(mdc.alpha));

    // Mu comes directly from the network as an unmolested value.
    // Note that this effectively means that the output layer of
    // the network should have an activation function at least as large as
    // the expected values.  It is best for the output
    // layer to be an IDENTITY activation function.
    //mdc.mu = mdc.mu;

    // Sigma comes from the network as an exponential in order to
    // ensure that it is positive-definite.
    mdc.sigma = Transforms.exp(mdc.sigma);

    return mdc;
}
 

/**
 * Local implementation of the E-step update of copy ratio posteriors
 *
 * @return a {@link SubroutineSignal} containing the update size (key: "error_norm")
 */
public SubroutineSignal updateCopyRatioPosteriorExpectationsLocal(final double admixingRatio) {
    /* step 1. fetch copy ratio emission data */
    final List<List<CoverageModelCopyRatioEmissionData>> copyRatioEmissionData = fetchCopyRatioEmissionDataLocal();

    /* step 2. run the forward-backward algorithm and calculate copy ratio posteriors */
    final INDArray sampleReadDepths = Transforms.exp(sampleMeanLogReadDepths, true);
    final List<CopyRatioExpectations> copyRatioPosteriorResults = sampleIndexStream()
            .parallel()
            .mapToObj(si -> copyRatioExpectationsCalculator.getCopyRatioPosteriorExpectations(
                    CopyRatioCallingMetadata.builder()
                            .sampleName(processedSampleNameList.get(si))
                            .sampleSexGenotypeData(processedSampleSexGenotypeData.get(si))
                            .sampleCoverageDepth(sampleReadDepths.getDouble(si))
                            .emissionCalculationStrategy(EmissionCalculationStrategy.HYBRID_POISSON_GAUSSIAN)
                            .build(),
                    processedTargetList,
                    copyRatioEmissionData.get(si)))
            .collect(Collectors.toList());

    /* update log chain posterior expectation */
    sampleLogChainPosteriors.assign(Nd4j.create(copyRatioPosteriorResults.stream()
            .mapToDouble(CopyRatioExpectations::getLogChainPosteriorProbability)
            .toArray(), new int[] {numSamples, 1}));

    /* sent the results back to workers */
    final ImmutablePair<INDArray, INDArray> copyRatioPosteriorDataPair =
            convertCopyRatioLatentPosteriorExpectationsToNDArray(copyRatioPosteriorResults);
    final INDArray log_c_st = copyRatioPosteriorDataPair.left;
    final INDArray var_log_c_st = copyRatioPosteriorDataPair.right;

    /* partition the pair of (log_c_st, var_log_c_st), sent the result to workers via broadcast-hash-map */
    pushToWorkers(mapINDArrayPairToBlocks(log_c_st.transpose(), var_log_c_st.transpose()),
            (p, cb) -> cb.cloneWithUpdatedCopyRatioPosteriors(
                    p.get(cb.getTargetSpaceBlock()).left.transpose(),
                    p.get(cb.getTargetSpaceBlock()).right.transpose(),
                    admixingRatio));
    cacheWorkers("after E-step update of copy ratio posteriors");

    /* collect subroutine signals */
    final List<SubroutineSignal> sigs = mapWorkersAndCollect(CoverageModelEMComputeBlock::getLatestMStepSignal);

    final double errorNormInfinity = Collections.max(sigs.stream()
            .map(sig -> sig.<Double>get(StandardSubroutineSignals.RESIDUAL_ERROR_NORM))
            .collect(Collectors.toList()));

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

private void calculate(INDArray labels, INDArray preOutput, IActivation activationFn, INDArray mask, INDArray scoreOutput, INDArray gradientOutput) {
    if (scoreOutput == null && gradientOutput == null) {
        throw new IllegalArgumentException("You have to provide at least one of scoreOutput or gradientOutput!");
    }
    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) + ") ");

    }
    labels = labels.castTo(preOutput.dataType());   //No-op if already correct dtype
    final INDArray postOutput = activationFn.getActivation(preOutput.dup(), true);

    final INDArray positive = labels;
    final INDArray negative = labels.eq(0.0).castTo(Nd4j.defaultFloatingPointType());
    final INDArray normFactor = negative.sum(true,1).castTo(Nd4j.defaultFloatingPointType()).muli(positive.sum(true,1));


    long examples = positive.size(0);
    for (int i = 0; i < examples; i++) {
        final INDArray locCfn = postOutput.getRow(i, true);
        final long[] shape = locCfn.shape();

        final INDArray locPositive = positive.getRow(i, true);
        final INDArray locNegative = negative.getRow(i, true);
        final Double locNormFactor = normFactor.getDouble(i);

        final int outSetSize = locNegative.sumNumber().intValue();
        if(outSetSize == 0 || outSetSize == locNegative.columns()){
            if (scoreOutput != null) {
                scoreOutput.getRow(i, true).assign(0);
            }

            if (gradientOutput != null) {
                gradientOutput.getRow(i, true).assign(0);
            }
        }else {
            final INDArray operandA = Nd4j.ones(shape[1], shape[0]).mmul(locCfn);
            final INDArray operandB = operandA.transpose();

            final INDArray pairwiseSub = Transforms.exp(operandA.sub(operandB));

            final INDArray selection = locPositive.transpose().mmul(locNegative);

            final INDArray classificationDifferences = pairwiseSub.muli(selection).divi(locNormFactor);

            if (scoreOutput != null) {
                if (mask != null) {
                    final INDArray perLabel = classificationDifferences.sum(0);
                    LossUtil.applyMask(perLabel, mask.getRow(i, true));
                    perLabel.sum(scoreOutput.getRow(i, true), 0);
                } else {
                    classificationDifferences.sum(scoreOutput.getRow(i, true), 0, 1);
                }
            }

            if (gradientOutput != null) {
                gradientOutput.getRow(i, true).assign(classificationDifferences.sum(true, 0).addi(classificationDifferences.sum(true,1).transposei().negi()));
            }
        }
    }

    if (gradientOutput != null) {
        gradientOutput.assign(activationFn.backprop(preOutput.dup(), gradientOutput).getFirst());
        //multiply with masks, always
        if (mask != null) {
            LossUtil.applyMask(gradientOutput, mask);
        }
    }
}
 

public MixtureDensityComponents extractComponents(INDArray output) {
    long outputSize = output.size(1);
    if (outputSize != (mLabelWidth + 2) * mMixtures) {
        throw new IllegalArgumentException(
                        "Network output size " + outputSize + " must be (labels+2)*mixtures where labels = "
                                        + mLabelWidth + " and mixtures = " + mMixtures);
    }

    MixtureDensityComponents mdc = new MixtureDensityComponents();

    // Output is 2 dimensional (samples, labels)
    //
    // For each label vector of length 'labels', we will have
    // an output vector of length '(labels + 2) * nMixtures.
    // The first nMixtures outputs will correspond to the 'alpha' for each mixture.
    // The second nMixtures outputs will correspond to the 'sigma' and the last nMixtures*labels
    // will correspond to the 'mu' (mean) of the output.

    // Reorganize these.
    // alpha = samples, 0 to nMixtures
    // mu = samples, nMixtures to 2*nMixtures
    // sigma = samples, 2*nMixtures to (labels + 2)*nMixtures
    // Alpha is then sub-divided through reshape by mixtures per label and samples.

    mdc.alpha = output.get(NDArrayIndex.all(), NDArrayIndex.interval(0, mMixtures));
    mdc.sigma = output.get(NDArrayIndex.all(), NDArrayIndex.interval(mMixtures, 2 * mMixtures));
    mdc.mu = output.get(NDArrayIndex.all(), NDArrayIndex.interval(2 * mMixtures, (mLabelWidth + 2) * mMixtures))
                    .reshape(output.size(0), mMixtures, mLabelWidth);

    // Alpha is a softmax because
    // the alpha should all sum to 1 for a given gaussian mixture.
    mdc.alpha = Nd4j.exec((CustomOp) new SoftMax(mdc.alpha, mdc.alpha, -1))[0];

    // Mu comes directly from the network as an unmolested value.
    // Note that this effectively means that the output layer of
    // the network should have an activation function at least as large as
    // the expected values.  It is best for the output
    // layer to be an IDENTITY activation function.
    //mdc.mu = mdc.mu;

    // Sigma comes from the network as an exponential in order to
    // ensure that it is positive-definite.
    mdc.sigma = Transforms.exp(mdc.sigma);

    return mdc;
}
 

public static INDArray activate(@NonNull INDArray boundingBoxPriors, @NonNull INDArray input, boolean nchw, LayerWorkspaceMgr layerWorkspaceMgr){
    if(!nchw)
        input = input.permute(0,3,1,2); //NHWC to NCHW

    long mb = input.size(0);
    long h = input.size(2);
    long w = input.size(3);
    long b = boundingBoxPriors.size(0);
    long c = input.size(1)/b-5;  //input.size(1) == b * (5 + C) -> C = (input.size(1)/b) - 5

    INDArray output = layerWorkspaceMgr.create(ArrayType.ACTIVATIONS, input.dataType(), input.shape(), 'c');
    INDArray output5 = output.reshape('c', mb, b, 5+c, h, w);
    INDArray output4 = output;  //output.get(all(), interval(0,5*b), all(), all());
    INDArray input4 = input.dup('c');    //input.get(all(), interval(0,5*b), all(), all()).dup('c');
    INDArray input5 = input4.reshape('c', mb, b, 5+c, h, w);

    //X/Y center in grid: sigmoid
    INDArray predictedXYCenterGrid = input5.get(all(), all(), interval(0,2), all(), all());
    Transforms.sigmoid(predictedXYCenterGrid, false);

    //width/height: prior * exp(input)
    INDArray predictedWHPreExp = input5.get(all(), all(), interval(2,4), all(), all());
    INDArray predictedWH = Transforms.exp(predictedWHPreExp, false);
    Broadcast.mul(predictedWH, boundingBoxPriors.castTo(input.dataType()), predictedWH, 1, 2);  //Box priors: [b, 2]; predictedWH: [mb, b, 2, h, w]

    //Confidence - sigmoid
    INDArray predictedConf = input5.get(all(), all(), point(4), all(), all());   //Shape: [mb, B, H, W]
    Transforms.sigmoid(predictedConf, false);

    output4.assign(input4);

    //Softmax
    //TODO OPTIMIZE?
    INDArray inputClassesPreSoftmax = input5.get(all(), all(), interval(5, 5+c), all(), all());   //Shape: [minibatch, C, H, W]
    INDArray classPredictionsPreSoftmax2d = inputClassesPreSoftmax.permute(0,1,3,4,2) //[minibatch, b, c, h, w] To [mb, b, h, w, c]
            .dup('c').reshape('c', new long[]{mb*b*h*w, c});
    Transforms.softmax(classPredictionsPreSoftmax2d, false);
    INDArray postSoftmax5d = classPredictionsPreSoftmax2d.reshape('c', mb, b, h, w, c ).permute(0, 1, 4, 2, 3);

    INDArray outputClasses = output5.get(all(), all(), interval(5, 5+c), all(), all());   //Shape: [minibatch, C, H, W]
    outputClasses.assign(postSoftmax5d);

    if(!nchw)
        output = output.permute(0,2,3,1);       //NCHW to NHWC

    return output;
}
 

@Override
public INDArray generateRandom(INDArray preOutDistributionParams) {
    INDArray gamma = activationFn.getActivation(preOutDistributionParams.dup(), false);

    INDArray lambda = Transforms.exp(gamma, true);

    //Inverse cumulative distribution function: -log(1-p)/lambda

    INDArray u = Nd4j.rand(preOutDistributionParams.shape());

    //Note here: if u ~ U(0,1) then 1-u ~ U(0,1)
    return Transforms.log(u, false).divi(lambda).negi();
}
 

/**
 * Calculate the reconstruction probability, as described in An & Cho, 2015 - "Variational Autoencoder based
 * Anomaly Detection using Reconstruction Probability" (Algorithm 4)<br>
 * The authors describe it as follows: "This is essentially the probability of the data being generated from a given
 * latent variable drawn from the approximate posterior distribution."<br>
 * <br>
 * Specifically, for each example x in the input, calculate p(x). Note however that p(x) is a stochastic (Monte-Carlo)
 * estimate of the true p(x), based on the specified number of samples. More samples will produce a more accurate
 * (lower variance) estimate of the true p(x) for the current model parameters.<br>
 * <br>
 * Internally uses {@link #reconstructionLogProbability(INDArray, int)} for the actual implementation.
 * That method may be more numerically stable in some cases.<br>
 * <br>
 * The returned array is a column vector of reconstruction probabilities, for each example. Thus, reconstruction probabilities
 * can (and should, for efficiency) be calculated in a batched manner.
 *
 * @param data       The data to calculate the reconstruction probability for
 * @param numSamples Number of samples with which to base the reconstruction probability on.
 * @return Column vector of reconstruction probabilities for each example (shape: [numExamples,1])
 */
public INDArray reconstructionProbability(INDArray data, int numSamples) {
    INDArray reconstructionLogProb = reconstructionLogProbability(data, numSamples);
    return Transforms.exp(reconstructionLogProb.castTo(DataType.DOUBLE), false);    //Cast to double to reduce risk of numerical underflow
}