def logp(self, x):
		n = self.n
		p = self.p
		s = self.s


		if s > 1:
			X = self._x_creation(x)

			result = self._normalizing_constant(n, p, s) + self._results_inner(n,x)
			return pm.dist_math.bound(result,
						tt.all(X <= 1), tt.all(X >= -1),
						self._check_pos_def(x),
						n > 0)
		else:
			X = x[self.tri_index]
			X = tt.fill_diagonal(X, 1)

			result = self._normalizing_constant(n, p, s)
			result += (n - 1.) * tt.log(tt.nlinalg.det(X))
			# n-1 probably needs to become structure[0]-1
			# I don't really know the likehood structure honestly

			return pm.dist_math.bound(result,
						tt.all(X <= 1), tt.all(X >= -1),
						matrix_pos_def(X),
						n > 0)

 def logp(self, x):
     n = self.n
     p = self.p
     # only defined for sum(p) == 1
     return bound(
         factln(n) + T.sum(x * T.log(p) - factln(x)), T.all(x >= 0), T.all(x <= n), T.eq(T.sum(x), n), n >= 0
     )

 def logp(self, x):
     n = self.n
     p = self.p
     # only defined for sum(p) == 1
     return bound(
         factln(n) + sum(x * log(p) - factln(x)),
         n >= 0,
         eq(sum(x), n),
         all(0 <= x), all(x <= n))

    def logp(self, value):
        k = self.k
        a = self.a

        # only defined for sum(value) == 1
        return bound(tt.sum(logpow(value, a - 1) - gammaln(a), axis=-1)
                     + gammaln(tt.sum(a, axis=-1)),
                     tt.all(value >= 0), tt.all(value <= 1),
                     k > 1, tt.all(a > 0))

    def logp(self, x):
        n = self.n
        p = self.p

        X = x[self.tri_index]
        X = t.fill_diagonal(X, 1)

        result = self._normalizing_constant(n, p)
        result += (n - 1.0) * log(det(X))
        return bound(result, n > 0, all(le(X, 1)), all(ge(X, -1)))

    def logp(self, value):
        n = self.n
        p = self.p

        return bound(factln(n) - factln(value).sum() + (value * tt.log(p)).sum(),
                     tt.all(value >= 0),
                     tt.all(0 <= p), tt.all(p <= 1),
                     tt.isclose(p.sum(), 1),
                     broadcast_conditions=False
        )

	def logp(self, x):
		n = self.n
		p = self.p
		s = self.s

		#X = x[self.tri_index] # need to correct
		#X = tt.fill_diagonal(X, 1) # need to correct
		X = self._x_creation(x)

		result = self._normalizing_constant(n, p, s) + self._results_inner(n,x)
		#result += (n - 1.) * T.log(det(X)) # n-1 probably needs to become structure[0]-1
		return pm.dist_math.bound(result,
					 tt.all(X <= 1), tt.all(X >= -1),
					 n > 0)

    def logp(self, x):
        n = self.n
        eta = self.eta

        X = x[self.tri_index]
        X = tt.fill_diagonal(X, 1)

        result = _lkj_normalizing_constant(eta, n)
        result += (eta - 1.) * tt.log(det(X))
        return bound(result,
                     tt.all(X <= 1), tt.all(X >= -1),
                     matrix_pos_def(X),
                     eta > 0,
                     broadcast_conditions=False
        )

    def __init__(self, a, transform=transforms.stick_breaking, *args, **kwargs):
        super(Dirichlet, self).__init__(transform=transform, *args, **kwargs)
        self.a = a
        self.k = a.shape[0]
        self.mean = a / sum(a)

        self.mode = switch(all(a > 1), (a - 1) / sum(a - 1), nan)

    def dlogp(inputs, gradients):
        g_logp, = gradients
        cov, delta = inputs

        g_logp.tag.test_value = floatX(1.)
        n, k = delta.shape

        chol_cov = cholesky(cov)
        diag = tt.nlinalg.diag(chol_cov)
        ok = tt.all(diag > 0)

        chol_cov = tt.switch(ok, chol_cov, tt.fill(chol_cov, 1))
        delta_trans = solve_lower(chol_cov, delta.T).T

        inner = n * tt.eye(k) - tt.dot(delta_trans.T, delta_trans)
        g_cov = solve_upper(chol_cov.T, inner)
        g_cov = solve_upper(chol_cov.T, g_cov.T)

        tau_delta = solve_upper(chol_cov.T, delta_trans.T)
        g_delta = tau_delta.T

        g_cov = tt.switch(ok, g_cov, -np.nan)
        g_delta = tt.switch(ok, g_delta, -np.nan)

        return [-0.5 * g_cov * g_logp, -g_delta * g_logp]

    def compute_weights(self, energies, attended_mask):
        """Compute weights from energies in softmax-like fashion.

        .. todo ::

            Use :class:`~blocks.bricks.Softmax`.

        Parameters
        ----------
        energies : :class:`~theano.Variable`
            The energies. Must be of the same shape as the mask.
        attended_mask : :class:`~theano.Variable`
            The mask for the attended. The index in the sequence must be
            the first dimension.

        Returns
        -------
        weights : :class:`~theano.Variable`
            Summing to 1 non-negative weights of the same shape
            as `energies`.

        """
        # Stabilize energies first and then exponentiate
        energies = energies - energies.max(axis=0)
        unnormalized_weights = tensor.exp(energies)
        if attended_mask:
            unnormalized_weights *= attended_mask

        # If mask consists of all zeros use 1 as the normalization coefficient
        normalization = (unnormalized_weights.sum(axis=0) +
                         tensor.all(1 - attended_mask, axis=0))
        return unnormalized_weights / normalization

    def logp(self, value):
        p_ = self.p
        k = self.k

        # Clip values before using them for indexing
        value_clip = tt.clip(value, 0, k - 1)

        p = p_ / tt.sum(p_, axis=-1, keepdims=True)

        if p.ndim > 1:
            pattern = (p.ndim - 1,) + tuple(range(p.ndim - 1))
            a = tt.log(p.dimshuffle(pattern)[value_clip])
        else:
            a = tt.log(p[value_clip])

        return bound(a, value >= 0, value <= (k - 1),
                     tt.all(p_ >= 0, axis=-1), tt.all(p <= 1, axis=-1))

    def in_transit(self, t, r=0.0, texp=None):
        """Get a list of timestamps that are in transit

        Args:
            t (vector): A vector of timestamps to be evaluated.
            r (Optional): The radii of the planets.
            texp (Optional[float]): The exposure time.

        Returns:
            The indices of the timestamps that are in transit.

        """

        z = tt.zeros_like(self.a)
        r = tt.as_tensor_variable(r) + z
        R = self.r_star + z

        # Wrap the times into time since transit
        hp = 0.5 * self.period
        dt = tt.mod(self._warp_times(t) - self.t0 + hp, self.period) - hp

        if self.ecc is None:
            # Equation 14 from Winn (2010)
            k = r / R
            arg = tt.square(1 + k) - tt.square(self.b)
            factor = R / (self.a * self.sin_incl)
            hdur = hp * tt.arcsin(factor * tt.sqrt(arg)) / np.pi
            t_start = -hdur
            t_end = hdur
            flag = z

        else:
            M_contact = self.contact_points_op(
                self.a, self.ecc, self.cos_omega, self.sin_omega,
                self.cos_incl + z, self.sin_incl + z, R + r)
            flag = M_contact[2]

            t_start = (M_contact[0] - self.M0) / self.n
            t_start = tt.mod(t_start + hp, self.period) - hp
            t_end = (M_contact[1] - self.M0) / self.n
            t_end = tt.mod(t_end + hp, self.period) - hp

            t_start = tt.switch(tt.gt(t_start, 0.0),
                                t_start - self.period, t_start)
            t_end = tt.switch(tt.lt(t_end, 0.0),
                              t_end + self.period, t_end)

        if texp is not None:
            t_start -= 0.5*texp
            t_end += 0.5*texp

        mask = tt.any(tt.and_(dt >= t_start, dt <= t_end), axis=-1)
        result = ifelse(tt.all(tt.eq(flag, 0)),
                        tt.arange(t.size)[mask],
                        tt.arange(t.size))

        return result

def negative_log_likelihood(actual, target):
    """
    :param actual: An (n_samples, n_labels) tensor where rows are normalized and actual[i,j] indicates the belief
        that on sample[i] the correct target is j.
    :param target: An (n_samples, ) tensor indicating the target label for each sample
    :return: The average (over samples) of the negative log-likelihood.
    """
    actual = tt.opt.assert_(actual, tt.all(abs(actual.sum(axis=1)-1) < 1e-7))  # Data must be normalized along axis 1.
    return negative_log_likelihood_dangerous(actual, target)

    def __init__(self, a, *args, **kwargs):
        super(Dirichlet, self).__init__(*args, **kwargs)
        self.a = a
        self.k = a.shape[0]
        self.mean = a / sum(a)

        self.mode = switch(all(a > 1),
                           (a - 1) / sum(a - 1),
                           nan)

    def logp(self, x):
        n = self.n
        p = self.p

        if x.ndim==2:
            x_sum = x.sum(axis=0)
            n_sum = n * x.shape[0]
        else:
            x_sum = x
            n_sum = n

        return bound(
            factln(n_sum) + tt.sum(x_sum * tt.log(p) - factln(x_sum)),
            tt.all(x >= 0),
            tt.all(x <= n),
            tt.eq(tt.sum(x_sum), n_sum),
            tt.isclose(p.sum(), 1),
            n >= 0)

    def __init__(self, a, transform=transforms.stick_breaking, *args, **kwargs):
        shape = a.shape[0]
        kwargs.setdefault("shape", shape)
        super(Dirichlet, self).__init__(transform=transform, *args, **kwargs)

        self.k = shape
        self.a = a
        self.mean = a / T.sum(a)

        self.mode = T.switch(T.all(a > 1), (a - 1) / T.sum(a - 1), np.nan)

	def logp(self, x):
		# x is assumed to be (s x n_elem) if s > 1 or n_elem
		n = self.n
		p = self.p
		s = self.s
		if s !=1:
			X = self._X_inner_creation(x)
			result = self._results_inner(n,p,s,x)
			return pm.dist_math.bound(result,
				tt.all(X <= 1), tt.all(X >= -1),
				n > 0)
		else:
			X = x[self.tri_index]
			X = tt.fill_diagonal(X, 1)
			result = self._normalizing_constant(n, p)
			result += (n - 1.) * tt.log(tt.nlinalg.det(X))
			return pm.dist_math.bound(result,
						 tt.all(X <= 1), tt.all(X >= -1),
						 n > 0)

    def __init__(self, a, transform=transforms.stick_breaking, *args, **kwargs):
        self.k = shape = a.shape[0]
        if "shape" not in kwargs.keys():
            kwargs.update({"shape": shape})
        super(Dirichlet, self).__init__(transform=transform, *args, **kwargs)
        self.a = a
        self.mean = a / sum(a)

        self.mode = switch(all(a > 1),
                           (a - 1) / sum(a - 1),
                           nan)

    def __init__(self, a, transform=transforms.stick_breaking,
                 *args, **kwargs):
        shape = a.shape[-1]
        kwargs.setdefault("shape", shape)
        super(Dirichlet, self).__init__(transform=transform, *args, **kwargs)

        self.k = tt.as_tensor_variable(shape)
        self.a = a = tt.as_tensor_variable(a)
        self.mean = a / tt.sum(a)

        self.mode = tt.switch(tt.all(a > 1),
                              (a - 1) / tt.sum(a - 1),
                              np.nan)