Python autograd.numpy.hstack() Examples

The following are code examples for showing how to use autograd.numpy.hstack(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: ISM_supervised_DR   Author: chieh-neu   File: kernel_lib.py    MIT License 6 votes vote down vote up
def Allocation_2_Y(allocation):
	N = np.size(allocation)
	unique_elements = np.unique(allocation)
	num_of_classes = len(unique_elements)
	class_ids = np.arange(num_of_classes)

	i = 0
	Y = np.zeros(num_of_classes)
	for m in allocation:
		class_label = np.where(unique_elements == m)[0]
		a_row = np.zeros(num_of_classes)
		a_row[class_label] = 1
		Y = np.hstack((Y, a_row))

	Y = np.reshape(Y, (N+1,num_of_classes))
	Y = np.delete(Y, 0, 0)

	return Y 
Example 2
Project: pilco-learner   Author: cryscan   File: gp.py    MIT License 6 votes vote down vote up
def log_pdf(self, hyp):
        x = np.atleast_2d(self.inputs)
        y = np.atleast_2d(self.targets)

        n, D = x.shape
        n, E = y.shape

        hyp = hyp.reshape(E, -1)
        K = self.kernel(hyp, x)  # [E, n, n]
        L = cholesky(K)
        alpha = np.hstack([solve(K[i], y[:, i]) for i in range(E)])
        y = y.flatten(order='F')

        logp = 0.5 * n * E * log(2 * np.pi) + 0.5 * np.dot(y, alpha) + np.sum(
            [log(np.diag(L[i])) for i in range(E)])

        return logp 
Example 3
Project: pilco-learner   Author: cryscan   File: control.py    MIT License 6 votes vote down vote up
def concat(con, sat, policy, m, s):
    max_u = policy.max_u
    E = len(max_u)
    D = len(m)

    F = D + E
    i, j = np.arange(D), np.arange(D, F)
    M = m
    S = fill_mat(s, np.zeros((F, F)))

    m, s, c = con(policy, m, s)
    M = np.hstack([M, m])
    S = fill_mat(s, S, j, j)
    q = np.matmul(S[np.ix_(i, i)], c)
    S = fill_mat(q, S, i, j)
    S = fill_mat(q.T, S, j, i)

    M, S, R = sat(M, S, j, max_u)
    C = np.hstack([np.eye(D), c]) @ R
    return M, S, C 
Example 4
Project: ceviche   Author: fancompute   File: fdfd.py    MIT License 6 votes vote down vote up
def _make_A(self, eps_vec):

        eps_vec_xx, eps_vec_yy = self._grid_average_2d(eps_vec)
        eps_vec_xx_inv = 1 / (eps_vec_xx + 1e-5)  # the 1e-5 is for numerical stability
        eps_vec_yy_inv = 1 / (eps_vec_yy + 1e-5)  # autograd throws 'divide by zero' errors.

        indices_diag = npa.vstack((npa.arange(self.N), npa.arange(self.N)))

        entries_DxEpsy,   indices_DxEpsy   = spsp_mult(self.entries_Dxb, self.indices_Dxb, eps_vec_yy_inv, indices_diag, self.N)
        entires_DxEpsyDx, indices_DxEpsyDx = spsp_mult(entries_DxEpsy, indices_DxEpsy, self.entries_Dxf, self.indices_Dxf, self.N)

        entries_DyEpsx,   indices_DyEpsx   = spsp_mult(self.entries_Dyb, self.indices_Dyb, eps_vec_xx_inv, indices_diag, self.N)
        entires_DyEpsxDy, indices_DyEpsxDy = spsp_mult(entries_DyEpsx, indices_DyEpsx, self.entries_Dyf, self.indices_Dyf, self.N)

        entries_d = 1 / EPSILON_0 * npa.hstack((entires_DxEpsyDx, entires_DyEpsxDy))
        indices_d = npa.hstack((indices_DxEpsyDx, indices_DyEpsxDy))

        entries_diag = MU_0 * self.omega**2 * npa.ones(self.N)

        entries_a = npa.hstack((entries_d, entries_diag))
        indices_a = npa.hstack((indices_d, indices_diag))

        return entries_a, indices_a 
Example 5
Project: ceviche   Author: fancompute   File: fdfd.py    MIT License 6 votes vote down vote up
def _solve_fn(self, eps_vec, entries_a, indices_a, Mz_vec):

        # convert the Mz current into Jx, Jy
        eps_vec_xx, eps_vec_yy = self._grid_average_2d(eps_vec)
        Jx_vec, Jy_vec = self._Hz_to_Ex_Ey(Mz_vec, eps_vec_xx, eps_vec_yy)

        # lump the current sources together and solve for electric field
        source_J_vec = npa.hstack((Jx_vec, Jy_vec))
        E_vec = sp_solve(entries_a, indices_a, source_J_vec)

        # strip out the x and y components of E and find the Hz component
        Ex_vec = E_vec[:self.N]
        Ey_vec = E_vec[self.N:]
        Hz_vec = self._Ex_Ey_to_Hz(Ex_vec, Ey_vec)

        return Ex_vec, Ey_vec, Hz_vec 
Example 6
Project: paragami   Author: rgiordan   File: pattern_containers.py    Apache License 2.0 6 votes vote down vote up
def flatten(self, folded_val, free=None, validate_value=None):
        free = self._free_with_default(free)
        valid, msg = self.validate_folded(
            folded_val, validate_value=validate_value)
        if not valid:
            raise ValueError(msg)

        # flat_length = self.flat_length(free)
        # offset = 0
        # flat_val = np.full(flat_length, float('nan'))
        flat_vals = []
        for pattern_name, pattern in self.__pattern_dict.items():
            pattern_flat_length = pattern.flat_length(free)
            # Containers must not mix free and non-free values, so do not
            # use default values for free.
            # flat_val[offset:(offset + pattern_flat_length)] = \
            flat_vals.append(
                pattern.flatten(
                    folded_val[pattern_name],
                    free=free,
                    validate_value=validate_value))
            #offset += pattern_flat_length
        return np.hstack(flat_vals) 
Example 7
Project: prediction-constrained-topic-models   Author: dtak   File: slda_utils__param_manager.py    MIT License 5 votes vote down vote up
def flatten_to_differentiable_param_vec(
        param_dict=None,
        topics_KV=None,
        w_CK=None,
        **unused_kwargs):
    """ Convert common parameters of sLDA into flat vector of reals.

    Examples
    --------
    >>> K = 2; V = 3; C = 2;
    >>> topics_KV = np.asarray([[0.6, 0.3, 0.1], [0.2, 0.1, 0.7]])
    >>> w_CK = np.asarray([[4.0, -4.0], [-1.0, 1.0]])
    >>> param_vec = flatten_to_differentiable_param_vec(
    ...     topics_KV=topics_KV,
    ...     w_CK=w_CK)
    >>> param_dict = unflatten_to_common_param_dict(
    ...     param_vec=param_vec, n_states=K, n_vocabs=V, n_labels=C)
    >>> print param_dict['w_CK']
    [[ 4. -4.]
     [-1.  1.]]

    >>> print param_dict['topics_KV']
    [[0.6 0.3 0.1]
     [0.2 0.1 0.7]]
    >>> np.allclose(param_dict['topics_KV'], topics_KV)
    True
    """
    if isinstance(param_dict, dict):
        topics_KV = param_dict['topics_KV']
        w_CK = param_dict['w_CK']
    return np.hstack([
        tfm__2D_rows_sum_to_one.to_diffable_arr(topics_KV).flatten(),
        w_CK.flatten()]) 
Example 8
Project: prediction-constrained-topic-models   Author: dtak   File: util_differentiable_transform__2D_rows_sum_to_one.py    MIT License 5 votes vote down vote up
def to_common_arr(
        log_topics_KVm1,
        min_eps=MIN_EPS,
        **kwargs):
    ''' Convert unconstrained topic weights to proper normalized topics

    Should handle any non-nan, non-inf input without numerical problems.

    Args
    ----
    log_topics_KVm1 : 2D array, size K x V-1

    Returns
    -------
    topics_KV : 2D array, size K x V
        minimum value of any entry will be min_eps
        each row will sum to 1.0 (+/- min_eps)
    '''
    K, Vm1 = log_topics_KVm1.shape
    V = Vm1 + 1
    log_topics_KV = np.hstack([
        log_topics_KVm1,
        np.zeros((K, 1))])
    log_topics_KV -= logsumexp(log_topics_KV, axis=1, keepdims=1)
    log_topics_KV += np.log1p(-V * min_eps)
    topics_KV = np.exp(log_topics_KV)

    return min_eps + topics_KV 
Example 9
Project: lightML   Author: jfzhang95   File: LinearClassification.py    MIT License 5 votes vote down vote up
def fit(self, X, y):
        # def loss function
        def calc_loss(W):
            y_pred = logistic_predictions(W, XMat)
            label_probabilities = y_pred * yMat + (1 - y_pred) * (1 - yMat)
            if self.reg is None:
                return -np.sum(np.log(label_probabilities))
            elif self.reg == 'l1':
                return -np.sum(np.log(label_probabilities))+np.sum(self.alpha*(np.abs(W[0:-1])))
            elif self.reg == 'l2':
                return -np.sum(np.log(label_probabilities))+np.sum(self.alpha*W[0:-1]*W[0:-1])
            else:
                print("the reg can only be None, l1 or l2!")
                return ValueError

        verbose = self.verbose
        print_step = self.print_step
        max_iters = self.max_iters

        XMat = np.array(X)
        yMat = np.array(y)

        if XMat.shape[0] != yMat.shape[0]:
            yMat = yMat.T
        assert XMat.shape[0] == yMat.shape[0]

        grad_fun = elementwise_grad(calc_loss)

        n_samples, n_features = X.shape
        n_outdim = y.shape[1]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])

        self.W = np.random.randn(n_features+1, n_outdim) * 0.1
        for it in range(max_iters + 1):
            loss = calc_loss(self.W)
            grad_params = grad_fun(self.W)
            # update params
            self.W -= self.lr * grad_params

            if verbose and it % print_step == 0:
                print('iteration %d / %d: loss %f' % (it, max_iters, loss)) 
Example 10
Project: lightML   Author: jfzhang95   File: LinearClassification.py    MIT License 5 votes vote down vote up
def predict(self, X):
        """Predict function"""
        XMat = np.array(X)
        n_samples = XMat.shape[0]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])
        ypred = logistic_predictions(self.W, XMat)
        for i in range(n_samples):
            if ypred[i] > 0.5:
                ypred[i] = 1
            else:
                ypred[i] = 0

        return ypred 
Example 11
Project: lightML   Author: jfzhang95   File: LinearRegression.py    MIT License 5 votes vote down vote up
def fit(self, X, y):
        # def loss function
        def calc_linear_loss(W):
            y_pred = np.dot(XMat, W)
            return np.sqrt((np.power(yMat - y_pred, 2))).mean()

        verbose = self.verbose
        print_step = self.print_step
        max_iters = self.max_iters

        XMat = np.array(X)
        yMat = np.array(y)

        if XMat.shape[0] != yMat.shape[0]:
            yMat = yMat.T
        assert XMat.shape[0] == yMat.shape[0]

        grad_fun = grad(calc_linear_loss)

        n_samples, n_features = X.shape
        n_outdim = y.shape[1]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])

        self.W = np.random.randn(n_features+1, n_outdim)
        for it in range(max_iters+1):
            loss = calc_linear_loss(self.W)
            grad_params = grad_fun(self.W)

            # update params
            self.W -= self.lr * grad_params

            if verbose and it % print_step == 0:
                print('iteration %d / %d: loss %f' % (it, max_iters, loss)) 
Example 12
Project: lightML   Author: jfzhang95   File: LinearRegression.py    MIT License 5 votes vote down vote up
def predict(self, X):
        """Predict function"""
        XMat = np.array(X)
        n_samples = XMat.shape[0]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])
        ypred = np.dot(XMat, self.W)

        return ypred 
Example 13
Project: lightML   Author: jfzhang95   File: LinearRegression.py    MIT License 5 votes vote down vote up
def fit(self, X, y):
        # def loss function
        def calc_linear_loss(W):
            y_pred = np.dot(XMat, W)
            return np.sqrt((np.power(yMat - y_pred, 2))).mean() \
                        + np.sum(self.alpha * W[0:-1] * W[0:-1])

        verbose = self.verbose
        print_step = self.print_step
        max_iters = self.max_iters

        XMat = np.array(X)
        yMat = np.array(y)

        if XMat.shape[0] != yMat.shape[0]:
            yMat = yMat.T
        assert XMat.shape[0] == yMat.shape[0]

        grad_fun = elementwise_grad(calc_linear_loss)

        n_samples, n_features = X.shape
        n_outdim = y.shape[1]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])

        self.W = np.random.randn(n_features+1, n_outdim)
        for it in range(max_iters+1):
            loss = calc_linear_loss(self.W)
            grad_params = grad_fun(self.W)

            # update params
            self.W -= self.lr * grad_params

            if verbose and it % print_step == 0:
                print('iteration %d / %d: loss %f' % (it, max_iters, loss)) 
Example 14
Project: lightML   Author: jfzhang95   File: LinearRegression.py    MIT License 5 votes vote down vote up
def fit(self, X, y):
        # def loss function
        def calc_linear_loss(W):
            y_pred = np.dot(XMat, W)
            return np.sqrt((np.power(yMat - y_pred, 2))).mean() \
                   + np.sum(self.alpha * np.abs(W[0:-1]))

        verbose = self.verbose
        print_step = self.print_step
        max_iters = self.max_iters

        XMat = np.array(X)
        yMat = np.array(y)

        if XMat.shape[0] != yMat.shape[0]:
            yMat = yMat.T
        assert XMat.shape[0] == yMat.shape[0]

        grad_fun = elementwise_grad(calc_linear_loss)

        n_samples, n_features = X.shape
        n_outdim = y.shape[1]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])

        self.W = np.random.randn(n_features + 1, n_outdim)
        for it in range(max_iters+1):
            loss = calc_linear_loss(self.W)
            grad_params = grad_fun(self.W)

            # update params
            self.W -= self.lr * grad_params

            if verbose and it % print_step == 0:
                print('iteration %d / %d: loss %f' % (it, max_iters, loss)) 
Example 15
Project: lightML   Author: jfzhang95   File: LinearRegression.py    MIT License 5 votes vote down vote up
def predict(self, X):
        """Predict function"""
        XMat = np.array(X)
        n_samples = XMat.shape[0]
        XMat = np.hstack([XMat, np.ones((n_samples, 1))])
        ypred = np.dot(XMat, self.W)
        return ypred 
Example 16
Project: autograd-forward   Author: BB-UCL   File: test_numpy.py    MIT License 5 votes vote down vote up
def test_hstack_1d(): combo_check(np.hstack, [0], [R(2), (R(2), R(2))]) 
Example 17
Project: autograd-forward   Author: BB-UCL   File: test_numpy.py    MIT License 5 votes vote down vote up
def test_hstack_2d(): combo_check(np.hstack, [0], [R(3, 2), (R(3, 4), R(3, 5))]) 
Example 18
Project: autograd-forward   Author: BB-UCL   File: test_numpy.py    MIT License 5 votes vote down vote up
def test_hstack_3d(): combo_check(np.hstack, [0], [R(2, 3, 4), (R(2, 1, 4), R(2, 5, 4))]) 
Example 19
Project: SyntheticStatistics   Author: BlissChapman   File: tst.py    MIT License 5 votes vote down vote up
def construct_z(X, Y, test_freqs, gaussian_width):
        """Construct the features Z to be used for testing with T^2 statistics.
        Z is defined in Eq.14 of Chwialkovski et al., 2015 (NIPS).

        test_freqs: J x d test frequencies

        Return a n x 2J numpy array. 2J because of sin and cos for each frequency.
        """
        if X.shape[0] != Y.shape[0]:
            raise ValueError('Sample size n must be the same for X and Y.')
        X = old_div(X,gaussian_width)
        Y = old_div(Y,gaussian_width)
        n, d = X.shape
        J = test_freqs.shape[0]
        # inverse Fourier transform (upto scaling) of the unit-width Gaussian kernel
        fx = np.exp(old_div(-np.sum(X**2, 1),2))[:, np.newaxis]
        fy = np.exp(old_div(-np.sum(Y**2, 1),2))[:, np.newaxis]
        # n x J
        x_freq = np.dot(X, test_freqs.T)
        y_freq = np.dot(Y, test_freqs.T)
        # zx: n x 2J
        zx = np.hstack((np.sin(x_freq)*fx, np.cos(x_freq)*fx))
        zy = np.hstack((np.sin(y_freq)*fy, np.cos(y_freq)*fy))
        z = zx-zy
        assert z.shape == (n, 2*J)
        return z 
Example 20
Project: SyntheticStatistics   Author: BlissChapman   File: data.py    MIT License 5 votes vote down vote up
def draw_mix_uniforms(pbs, pmix, n, seed):
        """
        Draw n samples from a mixture of uniform distributions whose bounds
        are specified in pbs, and mixing proportion is specified in pmix.

        Return a one-dimensional numpy array of n samples.
        """
        assert pbs.shape[0] == len(pmix)
        c = len(pmix)
        sam_list = []
        scale = pbs[:, 1] - pbs[:, 0]
        with util.NumpySeedContext(seed=seed):
            # counts for each mixture component
            counts = np.random.multinomial(n, pmix, size=1)
            # counts is a 2d array
            counts = counts[0]
            # For each component, draw from its corresponding uniform
            # distribution.
            for i, nc in enumerate(counts):
                #print i
                sam_i = stats.uniform.rvs(loc=pbs[i, 0], scale=scale[i], size=nc)
                #print 'pbs[{0},0]: {1}'.format(i, pbs[i, 0])
                #print 'sam_i: {0}'.format(sam_i)
                sam_list.append(sam_i)
            sample = np.hstack(sam_list)
            np.random.shuffle(sample)
        return sample

# end of SSMixUnif1D 
Example 21
Project: SyntheticStatistics   Author: BlissChapman   File: data.py    MIT License 5 votes vote down vote up
def sample(self, n, seed):
        rstate = np.random.get_state()
        np.random.seed(seed)

        d = self.d
        mean_y = np.hstack((self.my, np.zeros(d-1) ))
        X = np.random.randn(n, d)
        Y = np.random.randn(n, d) + mean_y
        np.random.set_state(rstate)
        return TSTData(X, Y, label='gmd_d%d'%self.d) 
Example 22
Project: ISM_supervised_DR   Author: chieh-neu   File: kernel_lib.py    MIT License 5 votes vote down vote up
def Y_2_allocation(Y):
	i = 0
	allocation = np.array([])
	for m in range(Y.shape[0]):
		allocation = np.hstack((allocation, np.where(Y[m] == 1)[0][0]))
		i += 1

	return allocation 
Example 23
Project: pilco-learner   Author: cryscan   File: cart_pole.py    MIT License 5 votes vote down vote up
def loss_cp(self, m, s):
    D0 = np.size(s, 1)
    D1 = D0 + 2 * len(self.angle)
    M = m
    S = s

    ell = self.p
    Q = np.dot(np.vstack([1, ell]), np.array([[1, ell]]))
    Q = fill_mat(Q, np.zeros((D1, D1)), [0, D0], [0, D0])
    Q = fill_mat(ell**2, Q, [D0 + 1], [D0 + 1])

    target = gaussian_trig(self.target, 0 * s, self.angle)[0]
    target = np.hstack([self.target, target])
    i = np.arange(D0)
    m, s, c = gaussian_trig(M, S, self.angle)
    q = np.dot(S[np.ix_(i, i)], c)
    M = np.hstack([M, m])
    S = np.vstack([np.hstack([S, q]), np.hstack([q.T, s])])

    w = self.width if hasattr(self, "width") else [1]
    L = np.array([0])
    S2 = np.array(0)
    for i in range(len(w)):
        self.z = target
        self.W = Q / w[i]**2
        r, s2, c = self.loss_sat(M, S)
        L = L + r
        S2 = S2 + s2

    return L / len(w) 
Example 24
Project: pilco-learner   Author: cryscan   File: util.py    MIT License 5 votes vote down vote up
def gaussian_trig(m, v, i, e=None):
    d = len(m)
    L = len(i)
    e = np.ones((1, L)) if e is None else np.atleast_2d(e)
    ee = np.vstack([e, e]).reshape(1, -1, order='F')

    mi = np.atleast_2d(m[i])
    vi = v[np.ix_(i, i)]
    vii = np.atleast_2d(np.diag(vi))

    M = np.vstack([e * exp(-vii / 2) * sin(mi), e * exp(-vii / 2) * cos(mi)])
    M = M.flatten(order='F')

    lq = -(vii.T + vii) / 2
    q = exp(lq)

    U1 = (exp(lq + vi) - q) * sin(mi.T - mi)
    U2 = (exp(lq - vi) - q) * sin(mi.T + mi)
    U3 = (exp(lq + vi) - q) * cos(mi.T - mi)
    U4 = (exp(lq - vi) - q) * cos(mi.T + mi)

    V = np.vstack(
        [np.hstack([U3 - U4, U1 + U2]),
         np.hstack([(U1 + U2).T, U3 + U4])])
    V = np.vstack([
        np.hstack([V[::2, ::2], V[::2, 1::2]]),
        np.hstack([V[1::2, ::2], V[1::2, 1::2]])
    ])
    V = np.dot(ee.T, ee) * V / 2

    C = np.hstack([np.diag(M[1::2]), -np.diag(M[::2])])
    C = np.hstack([C[:, ::2], C[:, 1::2]])
    C = fill_mat(C, np.zeros((d, 2 * L)), i, None)

    return M, V, C 
Example 25
Project: pilco-learner   Author: cryscan   File: util.py    MIT License 5 votes vote down vote up
def unwrap(p):
    return np.hstack([v.flatten() for v in p.values()]) 
Example 26
Project: pilco-learner   Author: cryscan   File: base.py    MIT License 5 votes vote down vote up
def train(gpmodel, plant, policy, x, y):
    Du = len(policy.max_u)
    dyni = np.asarray(plant.dyni)
    dyno = np.asarray(plant.dyno)
    difi = np.asarray(plant.difi)

    gpmodel.inputs = np.hstack([x[:, dyni], x[:, -Du:]])
    gpmodel.targets = y[:, dyno]
    gpmodel.targets[:, difi] = gpmodel.targets[:, difi] - x[:, dyno[difi]]
    gpmodel.optimize()

    hyp = gpmodel.hyp
    print(gpmodel.result['message'])
    print("Learned noise std:\n%s" % (str(np.exp(hyp[:, -1]))))
    print("SNRs:\n%s" % (str(np.exp(hyp[:, -2] - hyp[:, -1])))) 
Example 27
Project: ceviche   Author: fancompute   File: fdfd.py    MIT License 5 votes vote down vote up
def _make_A(self, eps_vec):

        C = - 1 / MU_0 * self.Dxf.dot(self.Dxb) \
            - 1 / MU_0 * self.Dyf.dot(self.Dyb)
        entries_c, indices_c = get_entries_indices(C)

        # indices into the diagonal of a sparse matrix
        entries_diag = - EPSILON_0 * self.omega**2 * eps_vec
        indices_diag = npa.vstack((npa.arange(self.N), npa.arange(self.N)))

        entries_a = npa.hstack((entries_diag, entries_c))
        indices_a = npa.hstack((indices_diag, indices_c))

        return entries_a, indices_a 
Example 28
Project: ceviche   Author: fancompute   File: fdfd.py    MIT License 5 votes vote down vote up
def _make_A(self, eps_vec):

        # notation: C = [[C11, C12], [C21, C22]]
        C11 = -1 / MU_0 * self.Dyf.dot(self.Dyb) 
        C22 = -1 / MU_0 * self.Dxf.dot(self.Dxb)
        C12 =  1 / MU_0 * self.Dyf.dot(self.Dxb)
        C21 =  1 / MU_0 * self.Dxf.dot(self.Dyb)

        # get entries and indices
        entries_c11, indices_c11 = get_entries_indices(C11)
        entries_c22, indices_c22 = get_entries_indices(C22)
        entries_c12, indices_c12 = get_entries_indices(C12)
        entries_c21, indices_c21 = get_entries_indices(C21)

        # shift the indices into each of the 4 quadrants
        indices_c22 += self.N       # shift into bottom right quadrant
        indices_c12[1,:] += self.N  # shift into top right quadrant
        indices_c21[0,:] += self.N  # shift into bottom left quadrant

        # get full matrix entries and indices
        entries_c = npa.hstack((entries_c11, entries_c12, entries_c21, entries_c22))
        indices_c = npa.hstack((indices_c11, indices_c12, indices_c21, indices_c22))

        # indices into the diagonal of a sparse matrix
        eps_vec_xx, eps_vec_yy, eps_vec_zz = self._grid_average_3d(eps_vec)
        entries_diag = - EPSILON_0 * self.omega**2 * npa.hstack((eps_vec_xx, eps_vec_yy))
        indices_diag = npa.vstack((npa.arange(2 * self.N), npa.arange(2 * self.N)))

        # put together the big A and return entries and indices
        entries_a = npa.hstack((entries_diag, entries_c))
        indices_a = npa.hstack((indices_diag, indices_c))
        return entries_a, indices_a 
Example 29
Project: kernel-gof   Author: wittawatj   File: ex2_prob_params.py    MIT License 5 votes vote down vote up
def job_mmd_opt(p, data_source, tr, te, r):
    """
    MMD test of Gretton et al., 2012 used as a goodness-of-fit test.
    Require the ability to sample from p i.e., the UnnormalizedDensity p has 
    to be able to return a non-None from get_datasource()

    With optimization. Gaussian kernel.
    """
    data = tr + te
    X = data.data()
    with util.ContextTimer() as t:
        # median heuristic 
        pds = p.get_datasource()
        datY = pds.sample(data.sample_size(), seed=r+294)
        Y = datY.data()
        XY = np.vstack((X, Y))

        med = util.meddistance(XY, subsample=1000)

        # Construct a list of kernels to try based on multiples of the median
        # heuristic
        #list_gwidth = np.hstack( (np.linspace(20, 40, 10), (med**2)
        #    *(2.0**np.linspace(-2, 2, 20) ) ) )
        list_gwidth = (med**2)*(2.0**np.linspace(-3, 3, 30) ) 
        list_gwidth.sort()
        candidate_kernels = [kernel.KGauss(gw2) for gw2 in list_gwidth]

        mmd_opt = mgof.QuadMMDGofOpt(p, n_permute=300, alpha=alpha, seed=r+56)
        mmd_result = mmd_opt.perform_test(data,
                candidate_kernels=candidate_kernels,
                tr_proportion=tr_proportion, reg=1e-3)
    return { 'test_result': mmd_result, 'time_secs': t.secs}


# Define our custom Job, which inherits from base class IndependentJob 
Example 30
Project: kernel-gof   Author: wittawatj   File: ex1_vary_n.py    MIT License 5 votes vote down vote up
def job_mmd_opt(p, data_source, tr, te, r):
    """
    MMD test of Gretton et al., 2012 used as a goodness-of-fit test.
    Require the ability to sample from p i.e., the UnnormalizedDensity p has 
    to be able to return a non-None from get_datasource()

    With optimization. Gaussian kernel.
    """
    data = tr + te
    X = data.data()
    with util.ContextTimer() as t:
        # median heuristic 
        pds = p.get_datasource()
        datY = pds.sample(data.sample_size(), seed=r+294)
        Y = datY.data()
        XY = np.vstack((X, Y))

        med = util.meddistance(XY, subsample=1000)

        # Construct a list of kernels to try based on multiples of the median
        # heuristic
        #list_gwidth = np.hstack( (np.linspace(20, 40, 10), (med**2)
        #    *(2.0**np.linspace(-2, 2, 20) ) ) )
        list_gwidth = (med**2)*(2.0**np.linspace(-4, 4, 30) ) 
        list_gwidth.sort()
        candidate_kernels = [kernel.KGauss(gw2) for gw2 in list_gwidth]

        mmd_opt = mgof.QuadMMDGofOpt(p, n_permute=300, alpha=alpha, seed=r)
        mmd_result = mmd_opt.perform_test(data,
                candidate_kernels=candidate_kernels,
                tr_proportion=tr_proportion, reg=1e-3)
    return { 'test_result': mmd_result, 'time_secs': t.secs} 
Example 31
Project: kernel-gof   Author: wittawatj   File: goftest.py    MIT License 5 votes vote down vote up
def optimize_auto_init(p, dat, J, **ops):
        """
        Optimize parameters by calling optimize_locs_widths(). Automatically 
        initialize the test locations and the Gaussian width.

        Return optimized locations, Gaussian width, optimization info
        """
        assert J>0
        # Use grid search to initialize the gwidth
        X = dat.data()
        n_gwidth_cand = 5
        gwidth_factors = 2.0**np.linspace(-3, 3, n_gwidth_cand) 
        med2 = util.meddistance(X, 1000)**2

        k = kernel.KGauss(med2*2)
        # fit a Gaussian to the data and draw to initialize V0
        V0 = util.fit_gaussian_draw(X, J, seed=829, reg=1e-6)
        list_gwidth = np.hstack( ( (med2)*gwidth_factors ) )
        besti, objs = GaussFSSD.grid_search_gwidth(p, dat, V0, list_gwidth)
        gwidth = list_gwidth[besti]
        assert util.is_real_num(gwidth), 'gwidth not real. Was %s'%str(gwidth)
        assert gwidth > 0, 'gwidth not positive. Was %.3g'%gwidth
        logging.info('After grid search, gwidth=%.3g'%gwidth)

        
        V_opt, gwidth_opt, info = GaussFSSD.optimize_locs_widths(p, dat,
                gwidth, V0, **ops) 

        # set the width bounds
        #fac_min = 5e-2
        #fac_max = 5e3
        #gwidth_lb = fac_min*med2
        #gwidth_ub = fac_max*med2
        #gwidth_opt = max(gwidth_lb, min(gwidth_opt, gwidth_ub))
        return V_opt, gwidth_opt, info 
Example 32
Project: kernel-gof   Author: wittawatj   File: goftest.py    MIT License 5 votes vote down vote up
def power_criterion(p, dat, b, c, test_locs, reg=1e-2):
        k = kernel.KIMQ(b=b, c=c)
        return FSSD.power_criterion(p, dat, k, test_locs, reg)

    #@staticmethod
    #def optimize_auto_init(p, dat, J, **ops):
    #    """
    #    Optimize parameters by calling optimize_locs_widths(). Automatically 
    #    initialize the test locations and the Gaussian width.

    #    Return optimized locations, Gaussian width, optimization info
    #    """
    #    assert J>0
    #    # Use grid search to initialize the gwidth
    #    X = dat.data()
    #    n_gwidth_cand = 5
    #    gwidth_factors = 2.0**np.linspace(-3, 3, n_gwidth_cand) 
    #    med2 = util.meddistance(X, 1000)**2

    #    k = kernel.KGauss(med2*2)
    #    # fit a Gaussian to the data and draw to initialize V0
    #    V0 = util.fit_gaussian_draw(X, J, seed=829, reg=1e-6)
    #    list_gwidth = np.hstack( ( (med2)*gwidth_factors ) )
    #    besti, objs = GaussFSSD.grid_search_gwidth(p, dat, V0, list_gwidth)
    #    gwidth = list_gwidth[besti]
    #    assert util.is_real_num(gwidth), 'gwidth not real. Was %s'%str(gwidth)
    #    assert gwidth > 0, 'gwidth not positive. Was %.3g'%gwidth
    #    logging.info('After grid search, gwidth=%.3g'%gwidth)

        
    #    V_opt, gwidth_opt, info = GaussFSSD.optimize_locs_widths(p, dat,
    #            gwidth, V0, **ops) 

    #    # set the width bounds
    #    #fac_min = 5e-2
    #    #fac_max = 5e3
    #    #gwidth_lb = fac_min*med2
    #    #gwidth_ub = fac_max*med2
    #    #gwidth_opt = max(gwidth_lb, min(gwidth_opt, gwidth_ub))
    #    return V_opt, gwidth_opt, info 
Example 33
Project: paragami   Author: rgiordan   File: simplex_patterns.py    Apache License 2.0 5 votes vote down vote up
def _unconstrain_simplex_jacobian(simplex_vec):
    """Get the unconstraining Jacobian for a single simplex vector.
    """
    return np.hstack(
        [ np.full(len(simplex_vec) - 1, -1 / simplex_vec[0])[:, None],
          np.diag(1 / simplex_vec[1: ]) ]) 
Example 34
Project: paragami   Author: rgiordan   File: simplex_patterns.py    Apache License 2.0 5 votes vote down vote up
def flat_indices(self, folded_bool, free=None):
        # If no indices are specified, save time and return an empty array.
        if not np.any(folded_bool):
            return np.array([], dtype=int)

        free = self._free_with_default(free)
        shape_ok, err_msg = self._validate_folded_shape(folded_bool)
        if not shape_ok:
            raise ValueError(err_msg)
        if not free:
            folded_indices = self.fold(
                np.arange(self.flat_length(False), dtype=int),
                validate_value=False, free=False)
            return folded_indices[folded_bool]
        else:
            # Every element of a particular simplex depends on all
            # the free values for that simplex.

            # The simplex is the last index, which moves the fastest.
            indices = []
            offset = 0
            free_simplex_length = self.__simplex_size - 1
            array_ranges = (range(n) for n in self.__array_shape)
            for ind in itertools.product(*array_ranges):
                if np.any(folded_bool[ind]):
                    free_inds = np.arange(
                        offset * free_simplex_length,
                        (offset + 1) * free_simplex_length,
                        dtype=int)
                    indices.append(free_inds)
                offset += 1
            if len(indices) > 0:
                return np.hstack(indices)
            else:
                return np.array([]) 
Example 35
Project: paragami   Author: rgiordan   File: pattern_containers.py    Apache License 2.0 5 votes vote down vote up
def flat_indices(self, folded_bool, free=None):
        free = self._free_with_default(free)
        valid, msg = self.validate_folded(folded_bool, validate_value=False)
        if not valid:
            raise ValueError(msg)

        flat_length = self.flat_length(free)
        offset = 0
        indices = []
        for pattern_name, pattern in self.__pattern_dict.items():
            pattern_flat_length = pattern.flat_length(free)
            # Containers must not mix free and non-free values, so do not
            # use default values for free.
            pattern_indices = pattern.flat_indices(
                folded_bool[pattern_name], free=free)
            if len(pattern_indices) > 0:
                indices.append(pattern_indices + offset)
            offset += pattern_flat_length
        if len(indices) > 0:
            return np.hstack(indices)
        else:
            return np.array([], dtype=int)


##########################
# An array of a pattern. 
Example 36
Project: paragami   Author: rgiordan   File: pattern_containers.py    Apache License 2.0 5 votes vote down vote up
def flatten(self, folded_val, free=None, validate_value=None):
        free = self._free_with_default(free)
        valid, msg = self.validate_folded(
            folded_val, validate_value=validate_value)
        if not valid:
            raise ValueError(msg)

        return np.hstack(np.array([
            self.__base_pattern.flatten(
                folded_val[item], free=free, validate_value=validate_value)
            for item in itertools.product(*self.__array_ranges)])) 
Example 37
Project: DistCal   Author: Srceh   File: GP_Beta_cal.py    Apache License 2.0 4 votes vote down vote up
def get_calibration(t_test, mu_test, sigma_test, mu_w, cov_w, mu_shift):

    n_y = numpy.shape(mu_test)[0]

    n_t = numpy.shape(t_test)[1]

    q_hat = numpy.zeros((n_y, n_t))

    s_hat = numpy.zeros((n_y, n_t))

    for i in range(0, n_y):

        ln_s = scipy.stats.norm.logpdf(x=t_test, loc=mu_test[i, :], scale=sigma_test[i, :]).reshape(-1, 1)

        feature_q = numpy.hstack([scipy.stats.norm.logcdf(x=t_test, loc=mu_test[i, :],
                                                          scale=sigma_test[i, :]).reshape(-1, 1),
                                  scipy.stats.norm.logsf(x=t_test, loc=mu_test[i, :],
                                                         scale=sigma_test[i, :]).reshape(-1, 1),
                                  numpy.ones((n_t, 1))])

        w_sample = scipy.stats.multivariate_normal.rvs(size=1024, mean=mu_w[i, :], cov=cov_w[i, :, :])

        w_sample[:, 0] = -numpy.exp(w_sample[:, 0] / mu_shift[0] + mu_shift[1])

        w_sample[:, 1] = numpy.exp(w_sample[:, 1] / mu_shift[2] + mu_shift[3])

        w_sample[:, 2] = w_sample[:, 2] / mu_shift[4] + mu_shift[5]

        raw_prod = numpy.matmul(feature_q, w_sample.transpose())

        MAX = raw_prod.copy()

        MAX[MAX < 0] = 0

        q_hat[i, :] = numpy.mean(numpy.exp(-MAX) / (numpy.exp(-MAX) + numpy.exp(raw_prod - MAX)), axis=1).ravel()

        tmp_de = numpy.where(raw_prod <= 0,
                             2 * numpy.log(1 + numpy.exp(raw_prod)),
                             2 * (raw_prod + numpy.log(1 + 1 / numpy.exp(raw_prod))))

        ln_s_hat = (raw_prod + numpy.log((w_sample[:, 0] + w_sample[:, 1]) * numpy.exp(feature_q[:, 0].reshape(-1, 1)) -
                                         w_sample[:, 0]) - feature_q[:, 0].reshape(-1, 1) -
                    feature_q[:, 1].reshape(-1, 1) - tmp_de) + ln_s

        mc_s_hat = numpy.exp(ln_s_hat)

        mc_s_hat[numpy.isnan(mc_s_hat)] = 0

        mc_s_hat[numpy.isinf(mc_s_hat)] = 0

        s_hat[i, :] = numpy.mean(mc_s_hat, axis=1).ravel()

    return s_hat, q_hat 
Example 38
Project: SyntheticStatistics   Author: BlissChapman   File: tst.py    MIT License 4 votes vote down vote up
def optimize_freqs_width(tst_data, alpha, n_test_freqs=10, max_iter=400,
            freqs_step_size=0.2, gwidth_step_size=0.01, batch_proportion=1.0,
            tol_fun=1e-3, seed=1):
        """Optimize the test frequencies and the Gaussian kernel width by
        maximizing the test power. X, Y should not be the same data as used
        in the actual test (i.e., should be a held-out set).

        - max_iter: #gradient descent iterations
        - batch_proportion: (0,1] value to be multipled with nx giving the batch
            size in stochastic gradient. 1 = full gradient ascent.
        - tol_fun: termination tolerance of the objective value

        Return (test_freqs, gaussian_width, info)
        """
        J = n_test_freqs
        """
        Optimize the empirical version of Lambda(T) i.e., the criterion used
        to optimize the test locations, for the test based
        on difference of mean embeddings with Gaussian kernel.
        Also optimize the Gaussian width.

        :return a theano function T |-> Lambda(T)
        """
        d = tst_data.dim()
        # set the seed
        rand_state = np.random.get_state()
        np.random.seed(seed)

        # draw frequencies randomly from the standard Gaussian.
        # TODO: Can we do better?
        T0 = np.random.randn(J, d)
        # reset the seed back to the original
        np.random.set_state(rand_state)

        # grid search to determine the initial gwidth
        mean_sd = tst_data.mean_std()
        scales = 2.0**np.linspace(-4, 3, 20)
        list_gwidth = np.hstack( (mean_sd*scales*(d**0.5), 2**np.linspace(-4, 4, 20) ))
        list_gwidth.sort()
        besti, powers = SmoothCFTest.grid_search_gwidth(tst_data, T0,
                list_gwidth, alpha)
        # initialize with the best width from the grid search
        gwidth0 = list_gwidth[besti]
        assert util.is_real_num(gwidth0), 'gwidth0 not real. Was %s'%str(gwidth0)
        assert gwidth0 > 0, 'gwidth0 not positive. Was %.3g'%gwidth0

        func_z = SmoothCFTest.construct_z_theano
        # info = optimization info
        T, gamma, info = optimize_T_gaussian_width(tst_data, T0, gwidth0, func_z,
                max_iter=max_iter, T_step_size=freqs_step_size,
                gwidth_step_size=gwidth_step_size, batch_proportion=batch_proportion,
                tol_fun=tol_fun)
        assert util.is_real_num(gamma), 'gamma is not real. Was %s' % str(gamma)

        ninfo = {'test_freqs': info['Ts'], 'test_freqs0': info['T0'],
                'gwidths': info['gwidths'], 'obj_values': info['obj_values'],
                'gwidth0': gwidth0, 'gwidth0_powers': powers}
        return (T, gamma, ninfo  ) 
Example 39
Project: SyntheticStatistics   Author: BlissChapman   File: tst.py    MIT License 4 votes vote down vote up
def optimize_locs_width(tst_data, alpha, n_test_locs=10, max_iter=400,
            locs_step_size=0.1, gwidth_step_size=0.01, batch_proportion=1.0,
            tol_fun=1e-3, reg=1e-5, seed=1):
        """Optimize the test locations and the Gaussian kernel width by
        maximizing the test power. X, Y should not be the same data as used
        in the actual test (i.e., should be a held-out set).

        - max_iter: #gradient descent iterations
        - batch_proportion: (0,1] value to be multipled with nx giving the batch
            size in stochastic gradient. 1 = full gradient ascent.
        - tol_fun: termination tolerance of the objective value

        Return (test_locs, gaussian_width, info)
        """
        J = n_test_locs
        """
        Optimize the empirical version of Lambda(T) i.e., the criterion used
        to optimize the test locations, for the test based
        on difference of mean embeddings with Gaussian kernel.
        Also optimize the Gaussian width.

        :return a theano function T |-> Lambda(T)
        """

        med = util.meddistance(tst_data.stack_xy(), 1000)
        T0 = MeanEmbeddingTest.init_locs_2randn(tst_data, n_test_locs,
                subsample=10000, seed=seed)
        #T0 = MeanEmbeddingTest.init_check_subset(tst_data, n_test_locs, med**2,
        #      n_cand=30, seed=seed+10)
        func_z = MeanEmbeddingTest.construct_z_theano
        # Use grid search to initialize the gwidth
        list_gwidth2 = np.hstack( ( (med**2) *(2.0**np.linspace(-3, 4, 30) ) ) )
        list_gwidth2.sort()
        besti, powers = MeanEmbeddingTest.grid_search_gwidth(tst_data, T0,
                list_gwidth2, alpha)
        gwidth0 = list_gwidth2[besti]
        assert util.is_real_num(gwidth0), 'gwidth0 not real. Was %s'%str(gwidth0)
        assert gwidth0 > 0, 'gwidth0 not positive. Was %.3g'%gwidth0

        # info = optimization info
        T, gamma, info = optimize_T_gaussian_width(tst_data, T0, gwidth0, func_z,
                max_iter=max_iter, T_step_size=locs_step_size,
                gwidth_step_size=gwidth_step_size, batch_proportion=batch_proportion,
                tol_fun=tol_fun, reg=reg)
        assert util.is_real_num(gamma), 'gamma is not real. Was %s' % str(gamma)

        ninfo = {'test_locs': info['Ts'], 'test_locs0': info['T0'],
                'gwidths': info['gwidths'], 'obj_values': info['obj_values'],
                'gwidth0': gwidth0, 'gwidth0_powers': powers}
        return (T, gamma, ninfo  ) 
Example 40
Project: networked-popularity   Author: avalanchesiqi   File: predictors.py    MIT License 4 votes vote down vote up
def arnet_predict(params, model_input, mode='train'):
    num_features, num_step = model_input.shape
    ar_coef = params[: -num_features + 1]
    num_input = len(ar_coef)
    link_weights = params[-num_features + 1:]

    yhat = None
    latent_interest = None

    if mode == 'train':
        for t in range(num_input, num_step):
            # AR part
            res = sum(ar_coef * model_input[0, t - num_input: t])
            if latent_interest is None:
                latent_interest = np.array([res])
            else:
                latent_interest = np.hstack((latent_interest, res))

            # network part
            for i in range(1, num_features):
                res += link_weights[i - 1] * model_input[i, t]

            # write in this way because auto grad does not support vector assignment
            if yhat is None:
                yhat = np.array([res])
            else:
                yhat = np.hstack((yhat, res))
    elif mode == 'test':
        input_views = model_input[0, :]
        for t in range(num_input):
            # AR part
            res = sum(ar_coef * input_views[-num_input:])
            if latent_interest is None:
                latent_interest = np.array([res])
            else:
                latent_interest = np.hstack((latent_interest, res))

            # network part
            for i in range(1, num_features):
                res += link_weights[i - 1] * model_input[i, t]

            # write in this way because auto grad does not support vector assignment
            if yhat is None:
                yhat = np.array([res])
            else:
                yhat = np.hstack((yhat, res))
            input_views = np.hstack((input_views, res))
    return yhat, latent_interest 
Example 41
Project: pilco-learner   Author: cryscan   File: base.py    MIT License 4 votes vote down vote up
def rollout(start, policy, plant, cost, H):
    """
    Generate a state trajectory using an ODE solver.
    """
    odei = plant.odei
    poli = plant.poli
    dyno = plant.dyno
    angi = plant.angi

    nX = len(odei)
    nU = len(policy.max_u)
    nA = len(angi)

    state = start
    x = np.zeros([H + 1, nX + 2 * nA])
    x[0, odei] = multivariate_normal(start, plant.noise)

    u = np.zeros([H, nU])
    y = np.zeros([H, nX])
    L = np.zeros(H)
    latent = np.zeros([H + 1, nX + nU])

    for i in range(H):
        s = x[i, odei]
        a, _, _ = gaussian_trig(s, 0 * np.eye(nX), angi)
        s = np.hstack([s, a])
        x[i, -2 * nA:] = s[-2 * nA:]

        if hasattr(policy, "fcn"):
            u[i, :], _, _ = policy.fcn(s[poli], 0 * np.eye(len(poli)))
        else:
            u[i, :] = policy.max_u * (2 * rand(nU) - 1)
        latent[i, :] = np.hstack([state, u[i, :]])

        dynamics = plant.dynamics
        dt = plant.dt
        next = odeint(dynamics, state[odei], [0, dt], args=(u[i, :], ))
        state = next[-1, :]
        x[i + 1, odei] = multivariate_normal(state[odei], plant.noise)

        if hasattr(cost, "fcn"):
            L[i] = cost.fcn(state[dyno], 0 * np.eye(len(dyno)))

    y = x[1:H + 1, :nX]
    x = np.hstack([x[:H, :], u[:H, :]])
    latent[H, :nX] = state

    return x, y, L, latent 
Example 42
Project: pilco-learner   Author: cryscan   File: gp.py    MIT License 4 votes vote down vote up
def gp0(self, m, s):
        """
        Compute joint predictions for MGP with uncertain inputs.
        """
        assert hasattr(self, "hyp")
        if not hasattr(self, "K"):
            self.cache()

        x = np.atleast_2d(self.inputs)
        y = np.atleast_2d(self.targets)
        n, D = x.shape
        n, E = y.shape

        X = self.hyp
        iK = self.iK
        beta = self.alpha

        m = np.atleast_2d(m)
        inp = x - m

        # Compute the predicted mean and IO covariance.
        iL = np.stack([np.diag(exp(-X[i, :D])) for i in range(E)])
        iN = np.matmul(inp, iL)
        B = iL @ s @ iL + np.eye(D)
        t = np.stack([solve(B[i].T, iN[i].T).T for i in range(E)])
        q = exp(-np.sum(iN * t, 2) / 2)
        qb = q * beta.T
        tiL = np.matmul(t, iL)
        c = exp(2 * X[:, D]) / sqrt(det(B))

        M = np.sum(qb, 1) * c
        V = (np.transpose(tiL, [0, 2, 1]) @ np.expand_dims(qb, 2)).reshape(
            E, D).T * c
        k = 2 * X[:, D].reshape(E, 1) - np.sum(iN**2, 2) / 2

        # Compute the predicted covariance.
        inp = np.expand_dims(inp, 0) / np.expand_dims(exp(2 * X[:, :D]), 1)
        ii = np.repeat(inp[:, newaxis, :, :], E, 1)
        ij = np.repeat(inp[newaxis, :, :, :], E, 0)

        iL = np.stack([np.diag(exp(-2 * X[i, :D])) for i in range(E)])
        siL = np.expand_dims(iL, 0) + np.expand_dims(iL, 1)
        R = np.matmul(s, siL) + np.eye(D)
        t = 1 / sqrt(det(R))
        iRs = np.stack(
            [solve(R.reshape(-1, D, D)[i], s) for i in range(E * E)])
        iRs = iRs.reshape(E, E, D, D)
        Q = exp(k[:, newaxis, :, newaxis] + k[newaxis, :, newaxis, :] +
                maha(ii, -ij, iRs / 2))

        S = np.einsum('ji,iljk,kl->il', beta, Q, beta)
        tr = np.hstack([np.sum(Q[i, i] * iK[i]) for i in range(E)])
        S = (S - np.diag(tr)) * t + np.diag(exp(2 * X[:, D]))
        S = S - np.matmul(M[:, newaxis], M[newaxis, :])

        return M, S, V 
Example 43
Project: kernel-gof   Author: wittawatj   File: ex3_vary_nlocs.py    MIT License 4 votes vote down vote up
def job_fssdq_opt(p, data_source, tr, te, r, J, null_sim=None):
    """
    FSSD with optimization on tr. Test on te. Use a Gaussian kernel.
    """
    if null_sim is None:
        null_sim = gof.FSSDH0SimCovObs(n_simulate=2000, seed=r)

    Xtr = tr.data()
    with util.ContextTimer() as t:
        # Use grid search to initialize the gwidth
        n_gwidth_cand = 5
        gwidth_factors = 2.0**np.linspace(-3, 3, n_gwidth_cand) 
        med2 = util.meddistance(Xtr, 1000)**2

        k = kernel.KGauss(med2*2)
        # fit a Gaussian to the data and draw to initialize V0
        V0 = util.fit_gaussian_draw(Xtr, J, seed=r+1, reg=1e-6)
        list_gwidth = np.hstack( ( (med2)*gwidth_factors ) )
        besti, objs = gof.GaussFSSD.grid_search_gwidth(p, tr, V0, list_gwidth)
        gwidth = list_gwidth[besti]
        assert util.is_real_num(gwidth), 'gwidth not real. Was %s'%str(gwidth)
        assert gwidth > 0, 'gwidth not positive. Was %.3g'%gwidth
        logging.info('After grid search, gwidth=%.3g'%gwidth)
        
        ops = {
            'reg': 1e-2,
            'max_iter': 50,
            'tol_fun': 1e-4,
            'disp': True,
            'locs_bounds_frac': 10.0,
            'gwidth_lb': 1e-1,
            'gwidth_ub': 1e3,
            }

        V_opt, gwidth_opt, info = gof.GaussFSSD.optimize_locs_widths(p, tr,
                gwidth, V0, **ops) 
        # Use the optimized parameters to construct a test
        k_opt = kernel.KGauss(gwidth_opt)
        fssd_opt = gof.FSSD(p, k_opt, V_opt, null_sim=null_sim, alpha=alpha)
        fssd_opt_result = fssd_opt.perform_test(te)
    return {'test_result': fssd_opt_result, 'time_secs': t.secs, 
            'goftest': fssd_opt, 'opt_info': info,
            } 
Example 44
Project: kernel-gof   Author: wittawatj   File: ex2_prob_params.py    MIT License 4 votes vote down vote up
def job_fssdJ1q_opt(p, data_source, tr, te, r, J=1, null_sim=None):
    """
    FSSD with optimization on tr. Test on te. Use a Gaussian kernel.
    """
    if null_sim is None:
        null_sim = gof.FSSDH0SimCovObs(n_simulate=2000, seed=r)

    Xtr = tr.data()
    with util.ContextTimer() as t:
        # Use grid search to initialize the gwidth
        n_gwidth_cand = 5
        gwidth_factors = 2.0**np.linspace(-3, 3, n_gwidth_cand) 
        med2 = util.meddistance(Xtr, 1000)**2

        k = kernel.KGauss(med2*2)
        # fit a Gaussian to the data and draw to initialize V0
        V0 = util.fit_gaussian_draw(Xtr, J, seed=r+1, reg=1e-6)
        list_gwidth = np.hstack( ( (med2)*gwidth_factors ) )
        besti, objs = gof.GaussFSSD.grid_search_gwidth(p, tr, V0, list_gwidth)
        gwidth = list_gwidth[besti]
        assert util.is_real_num(gwidth), 'gwidth not real. Was %s'%str(gwidth)
        assert gwidth > 0, 'gwidth not positive. Was %.3g'%gwidth
        logging.info('After grid search, gwidth=%.3g'%gwidth)
        
        ops = {
            'reg': 1e-2,
            'max_iter': 40,
            'tol_fun': 1e-4,
            'disp': True,
            'locs_bounds_frac':10.0,
            'gwidth_lb': 1e-1,
            'gwidth_ub': 1e4,
            }

        V_opt, gwidth_opt, info = gof.GaussFSSD.optimize_locs_widths(p, tr,
                gwidth, V0, **ops) 
        # Use the optimized parameters to construct a test
        k_opt = kernel.KGauss(gwidth_opt)
        fssd_opt = gof.FSSD(p, k_opt, V_opt, null_sim=null_sim, alpha=alpha)
        fssd_opt_result = fssd_opt.perform_test(te)
    return {'test_result': fssd_opt_result, 'time_secs': t.secs, 
            'goftest': fssd_opt, 'opt_info': info,
            } 
Example 45
Project: kernel-gof   Author: wittawatj   File: ex1_vary_n.py    MIT License 4 votes vote down vote up
def job_fssdJ1q_opt(p, data_source, tr, te, r, J=1, null_sim=None):
    """
    FSSD with optimization on tr. Test on te. Use a Gaussian kernel.
    """
    if null_sim is None:
        null_sim = gof.FSSDH0SimCovObs(n_simulate=2000, seed=r)

    Xtr = tr.data()
    with util.ContextTimer() as t:
        # Use grid search to initialize the gwidth
        n_gwidth_cand = 5
        gwidth_factors = 2.0**np.linspace(-3, 3, n_gwidth_cand) 
        med2 = util.meddistance(Xtr, 1000)**2

        k = kernel.KGauss(med2)
        # fit a Gaussian to the data and draw to initialize V0
        V0 = util.fit_gaussian_draw(Xtr, J, seed=r+1, reg=1e-6)
        list_gwidth = np.hstack( ( (med2)*gwidth_factors ) )
        besti, objs = gof.GaussFSSD.grid_search_gwidth(p, tr, V0, list_gwidth)
        gwidth = list_gwidth[besti]
        assert util.is_real_num(gwidth), 'gwidth not real. Was %s'%str(gwidth)
        assert gwidth > 0, 'gwidth not positive. Was %.3g'%gwidth
        logging.info('After grid search, gwidth=%.3g'%gwidth)
        
        ops = {
            'reg': 1e-2,
            'max_iter': 30,
            'tol_fun': 1e-5,
            'disp': True,
            'locs_bounds_frac':30.0,
            'gwidth_lb': 1e-1,
            'gwidth_ub': 1e4,
            }

        V_opt, gwidth_opt, info = gof.GaussFSSD.optimize_locs_widths(p, tr,
                gwidth, V0, **ops) 
        # Use the optimized parameters to construct a test
        k_opt = kernel.KGauss(gwidth_opt)
        fssd_opt = gof.FSSD(p, k_opt, V_opt, null_sim=null_sim, alpha=alpha)
        fssd_opt_result = fssd_opt.perform_test(te)
    return {'test_result': fssd_opt_result, 'time_secs': t.secs, 
            'goftest': fssd_opt, 'opt_info': info,
            } 
Example 46
Project: kernel-gof   Author: wittawatj   File: mmd.py    MIT License 4 votes vote down vote up
def perform_test(self, dat, candidate_kernels=None, return_mmdtest=False,
            tr_proportion=0.2, reg=1e-3):
        """
        dat: an instance of Data
        candidate_kernels: a list of Kernel's to choose from
        tr_proportion: proportion of sample to be used to choosing the best
            kernel
        reg: regularization parameter for the test power criterion 
        """
        with util.ContextTimer() as t:
            seed = self.seed
            p = self.p
            ds = p.get_datasource()
            p_sample = ds.sample(dat.sample_size(), seed=seed+77)
            xtr, xte = p_sample.split_tr_te(tr_proportion=tr_proportion, seed=seed+18)
            # ytr, yte are of type data.Data
            ytr, yte = dat.split_tr_te(tr_proportion=tr_proportion, seed=seed+12)

            # training and test data
            tr_tst_data = fdata.TSTData(xtr.data(), ytr.data())
            te_tst_data = fdata.TSTData(xte.data(), yte.data())

            if candidate_kernels is None:
                # Assume a Gaussian kernel. Construct a list of 
                # kernels to try based on multiples of the median heuristic
                med = util.meddistance(tr_tst_data.stack_xy(), 1000)
                list_gwidth = np.hstack( ( (med**2) *(2.0**np.linspace(-4, 4, 10) ) ) )
                list_gwidth.sort()
                candidate_kernels = [kernel.KGauss(gw2) for gw2 in list_gwidth]

            alpha = self.alpha

            # grid search to choose the best Gaussian width
            besti, powers = tst.QuadMMDTest.grid_search_kernel(tr_tst_data,
                    candidate_kernels, alpha, reg=reg)
            # perform test 
            best_ker = candidate_kernels[besti]
            mmdtest = tst.QuadMMDTest(best_ker, self.n_permute, alpha=alpha)
            results = mmdtest.perform_test(te_tst_data)
            if return_mmdtest:
                results['mmdtest'] = mmdtest

        results['time_secs'] = t.secs
        return results