def define_hparam_searchspace(infer_params):
    """Returns a list of hyperparameter sets. 

    Search space is defined by the following fields in :code:`infer_params`:

    - :code:`prior_indvdls`
    - :code:`cs`
    - :code:`L_cov_21s`
    - :code:`betas_indvdl`
    - :code:`betas_noise`

    Other fields are treated as constants.  

    All of the returned values, i.e. hyperparameters, are used to calculate 
    log marginal probabilities and the maximal value is taken. 

    :param InferParams infer_params: Parameters defining hyperparameter set.

    Returned value is a list of dicts, each of that is a set of 
    hyperparameters. It is created by the following code:

    .. code-block:: python

        hparamss = [
            {
                # Fixed hyperparameters
                'seed': seed, 
                'standardize': standardize, 
                'subtract_mu_reg': subtract_mu_reg, 
                'fix_mu_zero': fix_mu_zero, 
                'prior_var_mu': prior_var_mu, 
                'max_c': max_c, 
                'n_mc_samples': n_mc_samples, 
                'P_M1': P_M1, 
                'P_M2': P_M2, 
                'dist_noise': dist_noise, 
                'df_indvdl': df_indvdl, 
                'prior_scale': prior_scale, 

                # Varied hyperparameters
                'causality': causality, 
                'prior_indvdl': prior_indvdl, 
                'v_indvdl_1': v_indvdl_1, 
                'v_indvdl_2': v_indvdl_2, 
                'L_cov': L_cov, 
                'beta_coeff': beta_coeff, 
                'beta_noise': beta_noise
            }
            for prior_indvdl in prior_indvdls
            for v_indvdl_1 in cs
            for v_indvdl_2 in cs
            for L_cov in L_covs
            for causality in causalities
            for beta_coeff in (betas_indvdl if prior_indvdl == 'gg' else [None])
            for beta_noise in (betas_noise if dist_noise == 'gg' else [None])]

    :return: List of hyperparameter sets. 
    :rtype: list of dicts
    """
    assert(type(infer_params) == InferParams)
    return define_hparam_searchspace_main(
        seed=infer_params.seed, 
        standardize=infer_params.standardize, 
        subtract_mu_reg=infer_params.subtract_mu_reg, 
        fix_mu_zero=infer_params.fix_mu_zero, 
        prior_var_mu=infer_params.prior_var_mu, 
        max_c=infer_params.max_c, 
        n_mc_samples=infer_params.n_mc_samples, 
        P_M1=infer_params.P_M1, 
        P_M2=infer_params.P_M2, 
        df_indvdl=infer_params.df_indvdl, 
        dist_noise=infer_params.dist_noise, 
        prior_scale=infer_params.prior_scale, 
        prior_indvdls=infer_params.prior_indvdls, 
        cs=infer_params.cs, 
        L_cov_21s=infer_params.L_cov_21s, 
        betas_indvdl=infer_params.betas_indvdl, 
        betas_noise=infer_params.betas_noise, 
        causalities=infer_params.causalities, 
        scale_coeff=infer_params.scale_coeff
    )
 def ll_laplace(e):
    """Return of log likelihood of the standard Laplace distribution. 

    :param e: Sample values.
    :type e: ndarray, shape=(n_theta_sampling, n_samples)
    :return: pdf values at given samples. 
    :rtype: ndarray, shape=(n_theta_sampling, n_samples)
    """
    b = 1 / sqrt(2)
    return -log(2 * b) - abs(e) / b
 def ll_gg(e, beta):
    """Return of log likelihood of generalized Gaussian distributions. 
    """
    beta = float(beta)
    m = gamma(0.5 / beta) / ((2**(1 / beta)) * gamma(3 / (2 * beta)))
    
    return - 0.5 * power((e**2) / m, beta) + log(beta) \
           - gammaln(0.5 / beta) - (0.5 / beta) * log(2) - 0.5 * log(m)
 def sample_gg(scov, beta, n_samples, rng, dim, normalize):
    """Draw samples from GG distribution. 

    See: Gómez, E., Gomez-Viilegas, M. A., & Marin, J. M. (1998). 
    A multivariate generalization of the power exponential family of distributions. 
    Communications in Statistics-Theory and Methods, 27(3), 589-600.

    Parameters
    ----------
    scov: numpy.ndarray, shape=(dim, dim)
        A scaled covariance matrix. For beta=1, it becomes the covariance matrix. 

    beta: float
        Shape parameter, a positive value. 

    n_samples: int
        Number of samples. 

    rng: numpy.random.RandomState
        Random number generator. 

    dim: int 
        The dimension of the distribution. 

    normalize: bool
        If true, each dimension is normalized such that its variance is 1. 
    """
    # Draw samples from unit sphere
    gs = rng.normal(0, 1, size=(n_samples, dim))
    ss = gs + sign(sign(gs) + 1e-10) * 1e-10
    ns = np.linalg.norm(ss, ord=2, axis=1)[:, np.newaxis]
    us = ss / ns

    # Draw samples from GG
    S = np.linalg.cholesky(scov).T
    ts = rng.gamma(shape=dim / (2.0 * beta), scale=2.0, 
                   size=n_samples)[:, np.newaxis]
    xs = (ts**(1.0 / (2.0 * beta))) * us.dot(S)

    if normalize:
        v = (2**(1.0 / beta) * gamma((dim + 2.0) / (2.0 * beta))) / \
            (dim * gamma(dim / (2.0 * beta)))
        xs = xs / np.sqrt(v * np.diag(scov))

    return xs
 def multivariatet(mu, Sigma, N, M, rng):
    """Return a sample (or samples) from the multivariate t distribution.

    This function is adopted from 
    http://kennychowdhary.me/2013/03/python-code-to-generate-samples-from-multivariate-t/.

    :param mu: Mean.
    :type mu: ndarray, shape=(n_dim,), dtype=float
    :param Sigma: Scaling matrix.
    :type Sigma: ndarray, shape=(n_dim, n_dim), dtype=float
    :param float N: Degrees of freedom.
    :param int M: Number of samples to produce. 
    :param np.random.RandomState rng: Random number generator. 
    :return: M samples of (n_dim)-dimensional multivariate t distribution.
    :rtype: ndarray, shape=(n_samples, n_dim), dtype=float

    """
    d = len(Sigma)
    g = np.tile(rng.gamma(N/2., 2./N, M), (d, 1)).T
    Z = rng.multivariate_normal(np.zeros(d), Sigma, M)
    return mu + Z / np.sqrt(g)
 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
 def fetch_mu_indvdl(
    cache, prior_indvdl, L_cov_21, df_indvdl, beta_coeff):
    mu_indvdl1_ = cache[
        ('mu_indvdl1', prior_indvdl, L_cov_21, df_indvdl, beta_coeff)]
    mu_indvdl2_ = cache[
        ('mu_indvdl2', prior_indvdl, L_cov_21, df_indvdl, beta_coeff)]

    return mu_indvdl1_, mu_indvdl2_
 def fetch_mu1(cache, fix_mu_zero):
    return cache[('mu1', fix_mu_zero)]
 def fetch_mu2(cache, fix_mu_zero):
    return cache[('mu2', fix_mu_zero)]
 def fetch_h1(cache, prior_scale):
    return cache[('h1', prior_scale)]
 def fetch_h2(cache, prior_scale):
    return cache[('h2', prior_scale)]
 def fetch_b(cache):
    return cache['b']

from itertools import chain
from scipy.misc import logsumexp
from bmlingam.hparam import define_hparam_searchspace
from bmlingam.prob import ll_laplace, ll_gg, sample_gg, multivariatet
from bmlingam.utils import standardize_samples
from bmlingam.cache_mc import fetch_mu_indvdl, fetch_mu1, fetch_mu2, \
                              fetch_h1, fetch_h2, fetch_b
import numpy as np

        L = rng.multivariate_normal(
            np.zeros(L_cov.shape[0]), L_cov, n_mc_samples * n_samples)
    elif hparams['prior_indvdl'] == 'gg':
        beta = hparams['beta_coeff']
        L = sample_gg(L_cov, beta, n_mc_samples * n_samples, rng)
    else:
        raise ValueError('Invalid value of prior_indvdl: %s' %
            hparams['prior_indvdl'])
    del L_cov

    L_mu1 = L[:, 0].reshape(n_mc_samples, n_samples)
    L_mu2 = L[:, 1].reshape(n_mc_samples, n_samples)

    # Scale variance according to given hyperparameters and data std
    mu1_ = std_x[0] * hparams['v_indvdl_1'] * L_mu1
    mu2_ = std_x[1] * hparams['v_indvdl_2'] * L_mu2

    # ---- Sampling functions ----
    gen_randnones = lambda: (
        rng.randn(n_mc_samples, 1) * np.ones((1, n_samples)))
    gen_logn_ones = lambda: (
        rng.lognormal(size=(n_mc_samples, 1)) * np.ones((1, n_samples)))
    gen_uniform_ones = lambda: (
        rng.uniform(size=(n_mc_samples, 1)) * np.ones((1, n_samples)))

    # ---- Mean ----
    if hparams['fix_mu_zero']:
        mu1 = np.zeros((n_mc_samples, n_samples))
        mu2 = np.zeros((n_mc_samples, n_samples))
    elif hparams['prior_var_mu'] == 'auto':
        mu1 = np.sqrt(tau_cmmn[0]) * gen_randnones()
        mu2 = np.sqrt(tau_cmmn[1]) * gen_randnones()
    else:
        v = hparams['prior_var_mu']
        mu1 = np.sqrt(v) * gen_randnones()
        mu2 = np.sqrt(v) * gen_randnones()
        
    # ---- Noise variance ----
    if hparams['prior_scale'] == 'tr_normal':
        h1 = np.sqrt(tau_cmmn[0]) * np.abs(gen_randnones())
        h2 = np.sqrt(tau_cmmn[1]) * np.abs(gen_randnones())
    elif hparams['prior_scale'] == 'log_normal':
        h1 = np.sqrt(tau_cmmn[0]) * gen_logn_ones()
        h2 = np.sqrt(tau_cmmn[1]) * gen_logn_ones()
    elif hparams['prior_scale'] == 'uniform':
        h1 = np.sqrt(tau_cmmn[0]) * gen_uniform_ones()
        h2 = np.sqrt(tau_cmmn[1]) * gen_uniform_ones()
    else:
        raise ValueError(
            "Invalid value of prior_scale: %s" % hparams['prior_scale'])

    # ---- Regression coefficient ----
    b21 = np.sqrt(tau_cmmn[1]) * gen_randnones()
    b12 = np.sqrt(tau_cmmn[0]) * gen_randnones()
    del max_c, tau_cmmn

    # ---- Compute log probabilities of samples given parameters ----
    exp_nakami_xis_givenMgammatheta1 = \
        _comp_exp_nakami_xis_givenMgammatheta(
            1, n_mc_samples, xs, (mu1 + mu1_), (mu2 + mu2_), b21, None, h1, h2, 
            hparams)
    exp_nakami_xis_givenMgammatheta2 = \
        _comp_exp_nakami_xis_givenMgammatheta(
            2, n_mc_samples, xs, (mu1 + mu1_), (mu2 + mu2_), None, b12, h1, h2, 
            hparams)

    # ---- Averaging over MonteCarlo samples ----
    logps_M1 = (logsumexp(exp_nakami_xis_givenMgammatheta1, 0) -
                np.log(n_mc_samples))
    logps_M2 = (logsumexp(exp_nakami_xis_givenMgammatheta2, 0) - 
                np.log(n_mc_samples))
    logP_D_givenM1 = np.sum(logps_M1)    
    logP_D_givenM2 = np.sum(logps_M2)

    # ------ Log marginal probability (up to const) ------
    logP_M1_givenD = logP_D_givenM1 + np.log(hparams['P_M1'])
    logP_M2_givenD = logP_D_givenM2 + np.log(hparams['P_M2'])

    # ---- Sampling traces ----
    traces = {
        'mu1': mu1[:, 0], 
        'mu2': mu2[:, 0], 
        'h1': h1[:, 0], 
        'h2': h2[:, 0], 
        'mu1_': mu1_,
        'mu2_': mu2_, 
        'b12': b12[:, 0], 
        'b21': b21[:, 0], 
        'logps_M1': logps_M1, 
        'logps_M2': logps_M2
    }

    return logP_M1_givenD, logP_M2_givenD, traces

"""Supplemental functions
"""

def _comp_exp_nakami_xis_givenMgammatheta(
    M_id, n_theta_sampling, xs, mu1, mu2, b21, b12, h1, h2, hparams):
    """Compute probabilities of samples given parameters: 
    :math:`p(x_{l} - \\mu_{l} - b_{lm}x_{m}|{\\bf\\theta})`. 

    It is assumed a Laplace distribution as noise model. 

    :param int M_id: Model ID (1 or 2).
    :param int n_theta_sampling: Number of sampling iterations.
    :param xs: Observed samples. 
    :type xs: ndarray, shape=(n_samples, 2), dtype=float
    :param float mu1: Mean of the 1st signal.
    :param float mu2: Mean of the 2nd signal.
    :param float b21: Connection strength: 1 -> 2.
    :param float b12: Connection strength: 2 -> 1.
    :param float h1: Std of the 1st signal.
    :param float h2: Std of the 2nd signal.
    :return: Probabilities of samples given parameters. 
    :rtype: ndarray, shape=(n_theta_sampling, n_samples), dtype=float
    """
    gen_ones = lambda: np.ones((n_theta_sampling, 1))

    if hparams['dist_noise'] == 'laplace':

        logp = ll_laplace
