def _inverse_log_det_jacobian(self, y):
   y = self._maybe_assert_valid_y(y)
   event_dims = self._event_dims_tensor(y)
   return math_ops.reduce_sum(
       -math_ops.log1p(-y) +
       (1 / self.concentration - 1) * math_ops.log(-math_ops.log1p(-y)) +
       math_ops.log(self.scale / self.concentration),
       axis=event_dims)

 def _call_log_survival_function(self, value, name, **kwargs):
   with self._name_scope(name, values=[value]):
     value = ops.convert_to_tensor(value, name="value")
     try:
       return self._log_survival_function(value, **kwargs)
     except NotImplementedError:
       return math_ops.log1p(-self.cdf(value, **kwargs))

 def _forward(self, x):
   x = self._maybe_assert_valid_x(x)
   if self.power == 0.:
     return math_ops.exp(x)
   # If large x accuracy is an issue, consider using:
   # (1. + x * self.power)**(1. / self.power) when x >> 1.
   return math_ops.exp(math_ops.log1p(x * self.power) / self.power)

 def _inverse_log_det_jacobian(self, y):
   y = self._maybe_assert_valid(y)
   event_dims = self._event_dims_tensor(y)
   return math_ops.reduce_sum(
       math_ops.log(self.concentration1) + math_ops.log(self.concentration0) +
       (self.concentration1 - 1) * math_ops.log(y) +
       (self.concentration0 - 1) * math_ops.log1p(-y**self.concentration1),
       axis=event_dims)

 def _forward_log_det_jacobian(self, x):
   x = self._maybe_assert_valid_x(x)
   event_dims = self._event_dims_tensor(x)
   if self.power == 0.:
     return math_ops.reduce_sum(x, axis=event_dims)
   return (1. / self.power - 1.) * math_ops.reduce_sum(
       math_ops.log1p(x * self.power),
       axis=event_dims)

  def __init__(self,
               temperature,
               logits=None,
               probs=None,
               validate_args=False,
               allow_nan_stats=True,
               name="RelaxedBernoulli"):
    """Construct RelaxedBernoulli distributions.

    Args:
      temperature: An 0-D `Tensor`, representing the temperature
        of a set of RelaxedBernoulli distributions. The temperature should be
        positive.
      logits: An N-D `Tensor` representing the log-odds
        of a positive event. Each entry in the `Tensor` parametrizes
        an independent RelaxedBernoulli distribution where the probability of an
        event is sigmoid(logits). Only one of `logits` or `probs` should be
        passed in.
      probs: An N-D `Tensor` representing the probability of a positive event.
        Each entry in the `Tensor` parameterizes an independent Bernoulli
        distribution. Only one of `logits` or `probs` should be passed in.
      validate_args: Python `Boolean`, default `False`. When `True` distribution
        parameters are checked for validity despite possibly degrading runtime
        performance. When `False` invalid inputs may silently render incorrect
        outputs.
      allow_nan_stats: Python `Boolean`, default `True`. When `True`, statistics
        (e.g., mean, mode, variance) use the value "`NaN`" to indicate the
        result is undefined.  When `False`, an exception is raised if one or
        more of the statistic's batch members are undefined.
      name: `String` name prefixed to Ops created by this class.

    Raises:
      ValueError: If both `probs` and `logits` are passed, or if neither.
    """
    parameters = locals()
    with ops.name_scope(name, values=[logits, probs, temperature]) as ns:
      with ops.control_dependencies([check_ops.assert_positive(temperature)]
                                    if validate_args else []):
        self._temperature = array_ops.identity(temperature, name="temperature")

      self._logits, self._probs = distribution_util.get_logits_and_probs(
          logits=logits, probs=probs, validate_args=validate_args)
      dist = logistic.Logistic(self._logits / self._temperature,
                               1. / self._temperature,
                               validate_args=validate_args,
                               allow_nan_stats=allow_nan_stats,
                               name=ns)
      self._parameters = parameters

    def inverse_log_det_jacobian_fn(y):
      return -math_ops.log(y) - math_ops.log1p(-y)

    sigmoid_bijector = bijector.Inline(
        forward_fn=math_ops.sigmoid,
        inverse_fn=(lambda y: math_ops.log(y) - math_ops.log1p(-y)),
        inverse_log_det_jacobian_fn=inverse_log_det_jacobian_fn,
        name="sigmoid")
    super(RelaxedBernoulli, self).__init__(dist, sigmoid_bijector, name=name)

 def _forward_log_det_jacobian(self, x):
   # y = sinh((arcsinh(x) + skewness) * tailweight)
   # Using sinh' = cosh, arcsinh'(x) = 1 / sqrt(x**2 + 1),
   # dy/dx
   # = cosh((arcsinh(x) + skewness) * tailweight) * tailweight / sqrt(x**2 + 1)
   event_dims = self._event_dims_tensor(x)
   return math_ops.reduce_sum(
       log_cosh((arcsinh(x) + self.skewness) * self.tailweight) +
       math_ops.log(self.tailweight) - 0.5 * math_ops.log1p(x**2),
       axis=event_dims)

 def _inverse_log_det_jacobian(self, y):
   # x = sinh(arcsinh(y) / tailweight - skewness)
   # Using sinh' = cosh, arcsinh'(y) = 1 / sqrt(y**2 + 1),
   # dx/dy
   # = cosh(arcsinh(y) / tailweight - skewness)
   #     / (tailweight * sqrt(y**2 + 1))
   event_dims = self._event_dims_tensor(y)
   return math_ops.reduce_sum(
       log_cosh(arcsinh(y) / self.tailweight - self.skewness) -
       math_ops.log(self.tailweight) - 0.5 * math_ops.log1p(y**2),
       axis=event_dims)

 def _call_log_survival_function(self, value, name, **kwargs):
   with self._name_scope(name, values=[value]):
     value = _convert_to_tensor(
         value, name="value", preferred_dtype=self.dtype)
     try:
       return self._log_survival_function(value, **kwargs)
     except NotImplementedError as original_exception:
       try:
         return math_ops.log1p(-self.cdf(value, **kwargs))
       except NotImplementedError:
         raise original_exception

  def _log_prob(self, counts):
    if self.validate_args:
      counts = distribution_util.embed_check_nonnegative_discrete(
          counts, check_integer=True)
    counts *= array_ops.ones_like(self.probs)
    probs = self.probs * array_ops.ones_like(counts)

    safe_domain = array_ops.where(
        math_ops.equal(counts, 0.),
        array_ops.zeros_like(probs),
        probs)
    return counts * math_ops.log1p(-safe_domain) + math_ops.log(probs)

 def _cdf(self, x):
   if self.validate_args:
     x = distribution_util.embed_check_nonnegative_integer_form(x)
   else:
     # Whether or not x is integer-form, the following is well-defined.
     # However, scipy takes the floor, so we do too.
     x = math_ops.floor(x)
   x *= array_ops.ones_like(self.probs)
   return array_ops.where(
       x < 0.,
       array_ops.zeros_like(x),
       -math_ops.expm1((1. + x) * math_ops.log1p(-self.probs)))

 def _cdf(self, counts):
   if self.validate_args:
     # We set `check_integer=False` since the CDF is defined on whole real
     # line.
     counts = math_ops.floor(
         distribution_util.embed_check_nonnegative_discrete(
             counts, check_integer=False))
   counts *= array_ops.ones_like(self.probs)
   return array_ops.where(
       counts < 0.,
       array_ops.zeros_like(counts),
       -math_ops.expm1(
           (counts + 1) * math_ops.log1p(-self.probs)))

 def _log_prob(self, x):
   if self.validate_args:
     x = distribution_util.embed_check_nonnegative_integer_form(x)
   else:
     # For consistency with cdf, we take the floor.
     x = math_ops.floor(x)
   x *= array_ops.ones_like(self.probs)
   probs = self.probs * array_ops.ones_like(x)
   safe_domain = array_ops.where(
       math_ops.equal(x, 0.),
       array_ops.zeros_like(probs),
       probs)
   return x * math_ops.log1p(-safe_domain) + math_ops.log(probs)

 def _sample_n(self, n, seed=None):
   # Uniform variates must be sampled from the open-interval `(0, 1)` rather
   # than `[0, 1)`. To do so, we use `np.finfo(self.dtype.as_numpy_dtype).tiny`
   # because it is the smallest, positive, "normal" number. A "normal" number
   # is such that the mantissa has an implicit leading 1. Normal, positive
   # numbers x, y have the reasonable property that, `x + y >= max(x, y)`. In
   # this case, a subnormal number (i.e., np.nextafter) can cause us to sample
   # 0.
   uniform = random_ops.random_uniform(
       shape=array_ops.concat([[n], self.batch_shape_tensor()], 0),
       minval=np.finfo(self.dtype.as_numpy_dtype).tiny,
       maxval=1.,
       dtype=self.dtype,
       seed=seed)
   sampled = math_ops.log(uniform) - math_ops.log1p(-1. * uniform)
   return sampled * self.scale + self.loc

 def _sample_n(self, n, seed=None):
   shape = array_ops.concat([[n], self.batch_shape_tensor()], 0)
   # Uniform variates must be sampled from the open-interval `(-1, 1)` rather
   # than `[-1, 1)`. In the case of `(0, 1)` we'd use
   # `np.finfo(self.dtype.as_numpy_dtype).tiny` because it is the smallest,
   # positive, "normal" number. However, the concept of subnormality exists
   # only at zero; here we need the smallest usable number larger than -1,
   # i.e., `-1 + eps/2`.
   uniform_samples = random_ops.random_uniform(
       shape=shape,
       minval=np.nextafter(self.dtype.as_numpy_dtype(-1.),
                           self.dtype.as_numpy_dtype(0.)),
       maxval=1.,
       dtype=self.dtype,
       seed=seed)
   return (self.loc - self.scale * math_ops.sign(uniform_samples) *
           math_ops.log1p(-math_ops.abs(uniform_samples)))

def log_cdf_laplace(x, name="log_cdf_laplace"):
  """Log Laplace distribution function.

  This function calculates `Log[L(x)]`, where `L(x)` is the cumulative
  distribution function of the Laplace distribution, i.e.

  ```L(x) := 0.5 * int_{-infty}^x e^{-|t|} dt```

  For numerical accuracy, `L(x)` is computed in different ways depending on `x`,

  ```
  x <= 0:
    Log[L(x)] = Log[0.5] + x, which is exact

  0 < x:
    Log[L(x)] = Log[1 - 0.5 * e^{-x}], which is exact
  ```

  Args:
    x: `Tensor` of type `float32`, `float64`.
    name: Python string. A name for the operation (default="log_ndtr").

  Returns:
    `Tensor` with `dtype=x.dtype`.

  Raises:
    TypeError: if `x.dtype` is not handled.
  """

  with ops.name_scope(name, values=[x]):
    x = ops.convert_to_tensor(x, name="x")

    # For x < 0, L(x) = 0.5 * exp{x} exactly, so Log[L(x)] = log(0.5) + x.
    lower_solution = -np.log(2.) + x

    # safe_exp_neg_x = exp{-x} for x > 0, but is
    # bounded above by 1, which avoids
    #   log[1 - 1] = -inf for x = log(1/2), AND
    #   exp{-x} --> inf, for x << -1
    safe_exp_neg_x = math_ops.exp(-math_ops.abs(x))

    # log1p(z) = log(1 + z) approx z for |z| << 1. This approxmation is used
    # internally by log1p, rather than being done explicitly here.
    upper_solution = math_ops.log1p(-0.5 * safe_exp_neg_x)

    return array_ops.where(x < 0., lower_solution, upper_solution)

 def _log_unnormalized_prob(self, x):
   x = self._maybe_assert_valid_sample(x)
   return ((self.concentration1 - 1.) * math_ops.log(x)
           + (self.concentration0 - 1.) * math_ops.log1p(-x))

 def _inverse(self, y):
   y = self._maybe_assert_valid_y(y)
   return self.scale * (-math_ops.log1p(-y)) ** (1 / self.concentration)

def sigmoid_cross_entropy_with_logits(logits, targets, name=None):
  """Computes sigmoid cross entropy given `logits`.

  Measures the probability error in discrete classification tasks in which each
  class is independent and not mutually exclusive.  For instance, one could
  perform multilabel classification where a picture can contain both an elephant
  and a dog at the same time.

  For brevity, let `x = logits`, `z = targets`.  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)))

  `logits` and `targets` must have the same type and shape.

  Args:
    logits: A `Tensor` of type `float32` or `float64`.
    targets: A `Tensor` of the same type and shape as `logits`.
    name: A name for the operation (optional).

  Returns:
    A `Tensor` of the same shape as `logits` with the componentwise
    logistic losses.

  Raises:
    ValueError: If `logits` and `targets` do not have the same shape.
  """
  with ops.name_scope(name, "logistic_loss", [logits, targets]) as name:
    logits = ops.convert_to_tensor(logits, name="logits")
    targets = ops.convert_to_tensor(targets, name="targets")
    try:
      targets.get_shape().merge_with(logits.get_shape())
    except ValueError:
      raise ValueError("logits and targets must have the same shape (%s vs %s)"
                       % (logits.get_shape(), targets.get_shape()))

    # 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)))
    # To allow computing gradients at zero, we define custom versions of max and
    # abs functions.
    zeros = array_ops.zeros_like(logits, dtype=logits.dtype)
    cond = (logits >= zeros)
    relu_logits = array_ops.where(cond, logits, zeros)
    neg_abs_logits = array_ops.where(cond, -logits, logits)
    return math_ops.add(relu_logits - logits * targets,
                        math_ops.log1p(math_ops.exp(neg_abs_logits)),
                        name=name)

def weighted_cross_entropy_with_logits(targets, logits, pos_weight, name=None):
  """Computes a weighted cross entropy.

  This is like `sigmoid_cross_entropy_with_logits()` except that `pos_weight`,
  allows one to trade off recall and precision by up- or down-weighting the
  cost of a positive error relative to a negative error.

  The usual cross-entropy cost is defined as:

    targets * -log(sigmoid(logits)) + (1 - targets) * -log(1 - sigmoid(logits))

  The argument `pos_weight` is used as a multiplier for the positive targets:

    targets * -log(sigmoid(logits)) * pos_weight +
        (1 - targets) * -log(1 - sigmoid(logits))

  For brevity, let `x = logits`, `z = targets`, `q = pos_weight`.
  The loss is:

        qz * -log(sigmoid(x)) + (1 - z) * -log(1 - sigmoid(x))
      = qz * -log(1 / (1 + exp(-x))) + (1 - z) * -log(exp(-x) / (1 + exp(-x)))
      = qz * log(1 + exp(-x)) + (1 - z) * (-log(exp(-x)) + log(1 + exp(-x)))
      = qz * log(1 + exp(-x)) + (1 - z) * (x + log(1 + exp(-x))
      = (1 - z) * x + (qz +  1 - z) * log(1 + exp(-x))
      = (1 - z) * x + (1 + (q - 1) * z) * log(1 + exp(-x))

  Setting `l = (1 + (q - 1) * z)`, to ensure stability and avoid overflow,
  the implementation uses

      (1 - z) * x + l * (log(1 + exp(-abs(x))) + max(-x, 0))

  `logits` and `targets` must have the same type and shape.

  Args:
    targets: A `Tensor` of the same type and shape as `logits`.
    logits: A `Tensor` of type `float32` or `float64`.
    pos_weight: A coefficient to use on the positive examples.
    name: A name for the operation (optional).

  Returns:
    A `Tensor` of the same shape as `logits` with the componentwise
    weighted logistic losses.

  Raises:
    ValueError: If `logits` and `targets` do not have the same shape.
  """
  with ops.name_scope(name, "logistic_loss", [logits, targets]) as name:
    logits = ops.convert_to_tensor(logits, name="logits")
    targets = ops.convert_to_tensor(targets, name="targets")
    try:
      targets.get_shape().merge_with(logits.get_shape())
    except ValueError:
      raise ValueError("logits and targets must have the same shape (%s vs %s)"
                       % (logits.get_shape(), targets.get_shape()))

    # The logistic loss formula from above is
    #   (1 - z) * x + (1 + (q - 1) * z) * log(1 + exp(-x))
    # For x < 0, a more numerically stable formula is
    #   (1 - z) * x + (1 + (q - 1) * z) * log(1 + exp(x)) - l * x
    # To avoid branching, we use the combined version
    #   (1 - z) * x + l * (log(1 + exp(-abs(x))) + max(-x, 0))
    log_weight = 1 + (pos_weight - 1) * targets
    return math_ops.add(
        (1 - targets) * logits,
        log_weight * (math_ops.log1p(math_ops.exp(-math_ops.abs(logits))) +
                      nn_ops.relu(-logits)),
        name=name)