function sigmoidCrossEntropyWithLogits(labels, logits) {
	  if (!arraysEqual(labels.shape, logits.shape)) {
	    throw new ValueError("logits and labels must have the same shape, but got shapes " + (JSON.stringify(labels.shape) + " and " + JSON.stringify(logits.shape)));
	  }

	  return tidy(function () {
	    // The logistic loss formula from above is
	    //   x - x * z + log(1 + exp(-x))
	    // For x < 0, a more numerically stable formula is
	    //   -x * z + log(1 + exp(x))
	    // Note that these two expressions can be combined into the following:
	    //   max(x, 0) - x * z + log(1 + exp(-abs(x)))
	    var reluLogits = logits.relu();
	    var negAbsLogits = logits.abs().neg();
	    return reluLogits.sub(logits.mul(labels)).add(negAbsLogits.exp().log1p());
	  });
	}
From TensorFlow's implementation in nn_impl.py:

For brevity, let `x = logits`, `z = labels`.  The logistic loss is
     z * -log(sigmoid(x)) + (1 - z) * -log(1 - sigmoid(x))
   = z * -log(1 / (1 + exp(-x))) + (1 - z) * -log(exp(-x) / (1 + exp(-x)))
   = z * log(1 + exp(-x)) + (1 - z) * (-log(exp(-x)) + log(1 + exp(-x)))
   = z * log(1 + exp(-x)) + (1 - z) * (x + log(1 + exp(-x))
   = (1 - z) * x + log(1 + exp(-x))
   = x - x * z + log(1 + exp(-x))
For x < 0, to avoid overflow in exp(-x), we reformulate the above
     x - x * z + log(1 + exp(-x))
   = log(exp(x)) - x * z + log(1 + exp(-x))
   = - x * z + log(1 + exp(x))
Hence, to ensure stability and avoid overflow, the implementation uses this
equivalent formulation
   max(x, 0) - x * z + log(1 + exp(-abs(x)))

@param labels The labels.
@param logits The logits.