# Path: bmlingam/prob.py
# class GeneralizedGaussian(Continuous):
#     def __init__(self, mu=0.0, beta=None, cov=None, *args, **kwargs):
#         super(GeneralizedGaussian, self).__init__(*args, **kwargs)
#         # assert(mu.shape[0] == cov.shape[0] == cov.shape[1])
#         dim = mu.shape[0]
# 
#         self.mu = mu
#         self.beta = beta
#         self.prec = tt.nlinalg.pinv(cov)
#         # self.k = (dim * tt.gamma(dim / 2.0)) / \
#         #          ((np.pi**(dim / 2.0)) * tt.gamma(1 + dim / (2 * beta)) * (2**(1 + dim / (2 * beta))))
#         self.logk = tt.log(dim) + tt.gammaln(dim / 2.0) - \
#                     (dim / 2.0) * tt.log(np.pi) - \
#                     tt.gammaln(1 + dim / (2 * beta)) - \
#                     (1 + dim / (2 * beta)) * tt.log(2.0)
# 
#     def logp(self, value):
#         x = value - self.mu
# 
#         if x.tag.test_value.ndim == 1:
#             xpx = tt.sum(x * self.prec * x)
#             normalize = tt.log(self.prec)
#         else:
#             (x.dot(self.prec) * x).sum(axis=x.ndim - 1)
#             # normalize = tt.log(tt.nlinalg.Det(self.prec))
#             normalize = logabsdet(self.prec)
# 
#         # return tt.log(self.k) + 0.5 * normalize - 0.5 * xpx
#         return self.logk + 0.5 * normalize - 0.5 * xpx
# 
# Path: bmlingam/utils/standardize_samples.py
# def standardize_samples(xs, standardize):
#     xs = deepcopy(xs)
#     if standardize == 'keepratio':
#         s = np.std(xs)
#         xs = (xs - np.mean(xs, axis=0)) / s
#     elif standardize == 'scaling':
#         xs = xs / np.std(xs, axis=0)
#     elif standardize == 'commonscaling':
#         xs = xs / np.std(xs)
#     elif standardize is True:
#         xs = (xs - np.mean(xs, axis=0)) / np.std(xs, axis=0)
#     elif standardize is False:
#         xs = xs
#     else:
#         raise ValueError("Invalid value of standardize: %s" % standardize)
# 
#     return xs


import numpy as np
    import theano.tensor as tt
from parse import parse
from bmlingam.prob import GeneralizedGaussian
from bmlingam.utils import standardize_samples
    from pymc3 import Metropolis, sample
    from pymc3 import Normal, Laplace, StudentT, Model, HalfNormal, \
                      Deterministic, Gamma, Lognormal, Uniform
    from scipy.special import gamma
    from theano.tensor import sgn
    from theano import config
    from theano.tensor.slinalg import cholesky

def _get_L_cov(L_cov_21, floatX, Uniform, tt):
    if type(L_cov_21) in (float, int):
        return np.array([[1.0, L_cov_21], 
                         [L_cov_21, 1.0]]).astype(floatX)
    elif _is_uniform(L_cov_21):
        r = parse('U({:f},{:f})', L_cov_21.replace(' ', ''))
        L_cov_21_ = Uniform('L_cov_21_', lower=r[0], upper=r[1])
        return tt.stack([1.0, L_cov_21_, L_cov_21_, 1.0]).reshape((2, 2))

def _noise_model(hparams, h1, h2, xs_, Laplace, floatX, gamma):
    u"""Distribution of observation noise. 
    """
    # ---- Noise model ----
    if hparams['dist_noise'] == 'laplace':
        obs1 = lambda mu: Laplace(
            'x1s', mu=mu, b=h1 / np.float32(np.sqrt(2.)), 
            observed=xs_[:, 0], dtype=floatX
        )
        obs2 = lambda mu: Laplace(
            'x2s', mu=mu, b=h2 / np.float32(np.sqrt(2.)), 
            observed=xs_[:, 1], dtype=floatX
        )

    elif hparams['dist_noise'] == 'gg':
        beta = hparams['beta_noise']
        obs1 = lambda mu: GeneralizedGaussian(
            'x1s', mu=mu, beta=np.float32(beta), cov=np.eye(1), 
            observed=xs_[:, 0], dtype=floatX
        )

# Path: bmlingam/prob.py
# class GeneralizedGaussian(Continuous):
#     def __init__(self, mu=0.0, beta=None, cov=None, *args, **kwargs):
#         super(GeneralizedGaussian, self).__init__(*args, **kwargs)
#         # assert(mu.shape[0] == cov.shape[0] == cov.shape[1])
#         dim = mu.shape[0]
# 
#         self.mu = mu
#         self.beta = beta
#         self.prec = tt.nlinalg.pinv(cov)
#         # self.k = (dim * tt.gamma(dim / 2.0)) / \
#         #          ((np.pi**(dim / 2.0)) * tt.gamma(1 + dim / (2 * beta)) * (2**(1 + dim / (2 * beta))))
#         self.logk = tt.log(dim) + tt.gammaln(dim / 2.0) - \
#                     (dim / 2.0) * tt.log(np.pi) - \
#                     tt.gammaln(1 + dim / (2 * beta)) - \
#                     (1 + dim / (2 * beta)) * tt.log(2.0)
# 
#     def logp(self, value):
#         x = value - self.mu
# 
#         if x.tag.test_value.ndim == 1:
#             xpx = tt.sum(x * self.prec * x)
#             normalize = tt.log(self.prec)
#         else:
#             (x.dot(self.prec) * x).sum(axis=x.ndim - 1)
#             # normalize = tt.log(tt.nlinalg.Det(self.prec))
#             normalize = logabsdet(self.prec)
# 
#         # return tt.log(self.k) + 0.5 * normalize - 0.5 * xpx
#         return self.logk + 0.5 * normalize - 0.5 * xpx
# 
# Path: bmlingam/utils/standardize_samples.py
# def standardize_samples(xs, standardize):
#     xs = deepcopy(xs)
#     if standardize == 'keepratio':
#         s = np.std(xs)
#         xs = (xs - np.mean(xs, axis=0)) / s
#     elif standardize == 'scaling':
#         xs = xs / np.std(xs, axis=0)
#     elif standardize == 'commonscaling':
#         xs = xs / np.std(xs)
#     elif standardize is True:
#         xs = (xs - np.mean(xs, axis=0)) / np.std(xs, axis=0)
#     elif standardize is False:
#         xs = xs
#     else:
#         raise ValueError("Invalid value of standardize: %s" % standardize)
# 
#     return xs


# Path: bmlingam/bmlingam_pm3.py
import numpy as np
    import theano.tensor as tt
from parse import parse
from bmlingam.prob import GeneralizedGaussian
from bmlingam.utils import standardize_samples
    from pymc3 import Metropolis, sample
    from pymc3 import Normal, Laplace, StudentT, Model, HalfNormal, \
                      Deterministic, Gamma, Lognormal, Uniform
    from scipy.special import gamma
    from theano.tensor import sgn
    from theano import config
    from theano.tensor.slinalg import cholesky


def _get_L_cov(L_cov_21, floatX, Uniform, tt):
    if type(L_cov_21) in (float, int):
        return np.array([[1.0, L_cov_21], 
                         [L_cov_21, 1.0]]).astype(floatX)
    elif _is_uniform(L_cov_21):
        r = parse('U({:f},{:f})', L_cov_21.replace(' ', ''))
        L_cov_21_ = Uniform('L_cov_21_', lower=r[0], upper=r[1])
        return tt.stack([1.0, L_cov_21_, L_cov_21_, 1.0]).reshape((2, 2))

def _noise_model(hparams, h1, h2, xs_, Laplace, floatX, gamma):
    u"""Distribution of observation noise. 
    """
    # ---- Noise model ----
    if hparams['dist_noise'] == 'laplace':
        obs1 = lambda mu: Laplace(
            'x1s', mu=mu, b=h1 / np.float32(np.sqrt(2.)), 
            observed=xs_[:, 0], dtype=floatX
        )
        obs2 = lambda mu: Laplace(
            'x2s', mu=mu, b=h2 / np.float32(np.sqrt(2.)), 
            observed=xs_[:, 1], dtype=floatX
        )

    elif hparams['dist_noise'] == 'gg':
        beta = hparams['beta_noise']
        obs1 = lambda mu: GeneralizedGaussian(
            'x1s', mu=mu, beta=np.float32(beta), cov=np.eye(1), 
            observed=xs_[:, 0], dtype=floatX
        )

        obs2 = lambda mu: GeneralizedGaussian(