docs: public INDArray reconstructionLogProbability(INDArray data, int numSamples) {
          if (numSamples <= 0) {
              throw new IllegalArgumentException(
                              "Invalid input: numSamples must be > 0. Got: " + numSamples + " " + layerId());
          }
          if (reconstructionDistribution instanceof LossFunctionWrapper) {
              throw new UnsupportedOperationException("Cannot calculate reconstruction log probability when using "
                              + "a LossFunction (via LossFunctionWrapper) instead of a ReconstructionDistribution: ILossFunction "
                              + "instances are not in general probabilistic, hence it is not possible to calculate reconstruction probability "
                              + layerId());
          }
  
          data = data.castTo(dataType);
  
          //Forward pass through the encoder and mean for P(Z|X)
          LayerWorkspaceMgr workspaceMgr = LayerWorkspaceMgr.noWorkspaces();  //TODO add workspace support to this method
          setInput(data, workspaceMgr);
          VAEFwdHelper fwd = doForward(true, true, workspaceMgr);
          IActivation afn = layerConf().getActivationFn();
  
          //Forward pass through logStd^2 for P(Z|X)
          INDArray pzxLogStd2W = getParamWithNoise(VariationalAutoencoderParamInitializer.PZX_LOGSTD2_W, false, workspaceMgr);
          INDArray pzxLogStd2b = getParamWithNoise(VariationalAutoencoderParamInitializer.PZX_LOGSTD2_B, false, workspaceMgr);
  
          INDArray meanZ = fwd.pzxMeanPreOut;
          INDArray logStdev2Z = fwd.encoderActivations[fwd.encoderActivations.length - 1].mmul(pzxLogStd2W)
                          .addiRowVector(pzxLogStd2b);
          pzxActivationFn.getActivation(meanZ, false);
          pzxActivationFn.getActivation(logStdev2Z, false);
  
          INDArray pzxSigma = Transforms.exp(logStdev2Z, false);
          Transforms.sqrt(pzxSigma, false);
  
          val minibatch = input.size(0);
          val size = fwd.pzxMeanPreOut.size(1);
  
          INDArray pxzw = getParamWithNoise(VariationalAutoencoderParamInitializer.PXZ_W, false, workspaceMgr);
          INDArray pxzb = getParamWithNoise(VariationalAutoencoderParamInitializer.PXZ_B, false, workspaceMgr);
  
          INDArray[] decoderWeights = new INDArray[decoderLayerSizes.length];
          INDArray[] decoderBiases = new INDArray[decoderLayerSizes.length];
  
          for (int i = 0; i < decoderLayerSizes.length; i++) {
              String wKey = "d" + i + WEIGHT_KEY_SUFFIX;
              String bKey = "d" + i + BIAS_KEY_SUFFIX;
              decoderWeights[i] = getParamWithNoise(wKey, false, workspaceMgr);
              decoderBiases[i] = getParamWithNoise(bKey, false, workspaceMgr);
          }
  
          INDArray sumReconstructionNegLogProbability = null;
          for (int i = 0; i < numSamples; i++) {
              INDArray e = Nd4j.randn(dataType, minibatch, size);
              INDArray z = e.muli(pzxSigma).addi(meanZ); //z = mu + sigma * e, with e ~ N(0,1)
  
              //Do forward pass through decoder
              int nDecoderLayers = decoderLayerSizes.length;
              INDArray currentActivations = z;
              for (int j = 0; j < nDecoderLayers; j++) {
                  currentActivations = currentActivations.mmul(decoderWeights[j]).addiRowVector(decoderBiases[j]);
                  afn.getActivation(currentActivations, false);
              }
  
              //And calculate reconstruction distribution preOut
              INDArray pxzDistributionPreOut = currentActivations.mmul(pxzw).addiRowVector(pxzb);
  
              if (i == 0) {
                  sumReconstructionNegLogProbability =
                                  reconstructionDistribution.exampleNegLogProbability(data, pxzDistributionPreOut);
              } else {
                  sumReconstructionNegLogProbability
                                  .addi(reconstructionDistribution.exampleNegLogProbability(data, pxzDistributionPreOut));
              }
          }
  
          setInput(null, workspaceMgr);
          return sumReconstructionNegLogProbability.divi(-numSamples);
      }