Gradient of diagonal of covariance with respect to parameters.
973279
def dKdiag_dtheta(self, dL_dKdiag, X, target):
        """"""
        target += 2.*self.mapping.df_dtheta(dL_dKdiag[:, None]*self.mapping.f(X), X)
1095379
def _covar_mstep_spherical(*args):
    """"""
    cv = _covar_mstep_diag(*args)
    return np.tile(cv.mean(axis=1)[:, np.newaxis], (1, cv.shape[1]))
453888
def gradient(self):
        """"""
        return (self.functional.gradient +
                (2 * self.quadratic_coeff) * IdentityOperator(self.domain) +
                ConstantOperator(self.linear_term))
1327593
def _msge_with_gradient_overdetermined(data, delta, xvschema, skipstep, p):
    
    t, m, l = data.shape
    d = None
    l, k = 0, 0
    nt = np.ceil(t / skipstep)
    for trainset, testset in xvschema(t, skipstep):

        a, b = _construct_var_eqns(atleast_3d(data[trainset, :, :]), p)
        c, d = _construct_var_eqns(atleast_3d(data[testset, :, :]), p)

        e = sp.linalg.inv(np.eye(a.shape[1]) * delta ** 2 + a.T.dot(a))

        ba = b.transpose().dot(a)
        dc = d.transpose().dot(c)
        bae = ba.dot(e)
        baee = bae.dot(e)
        baecc = bae.dot(c.transpose().dot(c))

        l += np.sum(baecc * bae - 2 * bae * dc) + np.sum(d ** 2)
        k += np.sum(baee * dc - baecc * baee) * 4 * delta

    return l / (nt * d.size), k / (nt * d.size)
1361461
def eighCovarianceMatrix(cfile):
    
    # precompute eigenvalue decomposition
    K = np.loadtxt(cfile+'.cov')
    K+= 1e-4*sp.eye(K.shape[0])
    S,U = la.eigh(K); S=S[::-1]; U=U[:,::-1]
    np.savetxt(cfile+'.cov.eval',S,fmt='%.6f')
    np.savetxt(cfile+'.cov.evec',U,fmt='%.6f')
1146495
def normalizedGraylevelVariance(img):
    '''''''
    mean, stdev = cv2.meanStdDev(img)
    s = stdev[0]**2 / mean[0]
    return s[0]
675628
def gradients_X_X2(self, dL_dK, X, X2):
        """"""
        return self._comp_grads(dL_dK, X, X2)[3:]
1195563
def fgrad_y_psi(self, y, return_covar_chain=False):
        
        mpsi = self.psi

        w, s, r, d = self.fgrad_y(y, return_precalc=True)
        gradients = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
        for i in range(len(mpsi)):
            a,b,c  = mpsi[i]
            gradients[:, :, i, 0] = (b * (1.0/np.cosh(s[i])) ** 2).T
            gradients[:, :, i, 1] = a * (d[i] - 2.0 * s[i] * r[i] * (1.0/np.cosh(s[i])) ** 2).T
            gradients[:, :, i, 2] = (-2.0 * a * (b ** 2) * r[i] * ((1.0 / np.cosh(s[i])) ** 2)).T
        gradients[:, :, 0, 3] = 1.0

        if return_covar_chain:
            covar_grad_chain = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
            for i in range(len(mpsi)):
                a,b,c = mpsi[i]
                covar_grad_chain[:, :, i, 0] = (r[i]).T
                covar_grad_chain[:, :, i, 1] = (a * (y + c) * ((1.0 / np.cosh(s[i])) ** 2).T)
                covar_grad_chain[:, :, i, 2] = a * b * ((1.0 / np.cosh(s[i])) ** 2).T
            covar_grad_chain[:, :, 0, 3] = y
            return gradients, covar_grad_chain

        return gradients
515222
def _forward_step(self, Y, T_X, Var_X, Var_dX, rho_e, sigma2_e, weight):
        
        # We currently only implement diagonal form
        # of covariance matrix for Var_X, Var_dX and T_X, which means
        # each dimension of X is independent and their innovation noise
        # are also independent. Note that log_p_data takes this assumption.
        if Var_X.ndim == 1:
            inv_Var_X = np.diag(1 / Var_X)
            half_log_det_Var_X = np.sum(np.log(Var_X)) / 2.0
            Var_X = np.diag(Var_X)
            # the marginal variance of X
        else:
            half_log_det_Var_X = self._half_log_det(Var_X)
            inv_Var_X = np.linalg.inv(Var_X)
        if Var_dX.ndim == 1:
            inv_Var_dX = np.diag(1 / Var_dX)
            half_log_det_Var_dX = np.sum(np.log(Var_dX)) / 2.0
            Var_dX = np.diag(Var_dX)
            # the marginal variance of Delta X (the change of X from
            # previous time point)
        else:
            inv_Var_dX = np.linalg.inv(Var_dX)
            half_log_det_Var_dX = self._half_log_det(Var_dX)
        if T_X.ndim == 1:
            T_X = np.diag(T_X)
            # Transfer function of X: the expected mean of X at t+1
            # time point is T_x * X
        [n_T, n_V] = np.shape(Y)
        # numbers of time points and voxels
        mu = [None] * n_T
        # posterior mean of X, conditioned on all data up till the current
        # time point
        Gamma_inv = [None] * n_T
        # inverse of poterior Gamma.
        mu_Gamma_inv = [None] * n_T
        # mu * inv(Gamma)
        log_p_data = - np.log(np.pi * 2) * (n_T * n_V) / 2 \
            - half_log_det_Var_X - np.sum(np.log(sigma2_e)) * n_T / 2.0\
            + np.sum(np.log(1 - rho_e**2)) / 2.0 - half_log_det_Var_dX \
            * (n_T - 1)
        # This is the term to be incremented by c_n at each time step.
        # We first add all the fixed terms to it.

        # The following are a few fixed terms.
        Lambda_0 = np.dot(T_X, np.dot(inv_Var_dX, T_X.T)) \
            + np.dot(weight * rho_e**2 / sigma2_e, weight.T)
        H = np.dot(inv_Var_dX, T_X.T) + np.dot(weight * rho_e / sigma2_e,
                                               weight.T)
        Lambda_1 = inv_Var_dX + np.dot(weight / sigma2_e, weight.T)

        Gamma_inv[0] = inv_Var_X + np.dot(
            weight * (1 - rho_e**2) / sigma2_e, weight.T)
        # We might not need this and only use linalg.solve for related terms.
        mu_Gamma_inv[0] = np.dot(
            Y[0, :] * (1 - rho_e**2) / sigma2_e, weight.T)
        mu[0] = np.linalg.solve(Gamma_inv[0], mu_Gamma_inv[0])
        log_p_data -= 0.5 * np.sum(Y[0, :]**2 * (1 - rho_e**2) / sigma2_e)
        # This is the term added for the first time point.

        deltaY = Y[1:, :] - rho_e * Y[:-1, :]
        deltaY_sigma2inv_rho_weightT = np.dot(
            deltaY / sigma2_e * rho_e, weight.T)
        for t in np.arange(1, n_T):
            Gamma_tilde_inv = Lambda_0 + Gamma_inv[t - 1]
            tmp = np.linalg.solve(Gamma_tilde_inv, H.T)
            Gamma_inv[t] = Lambda_1 - np.dot(H, tmp)
            mu_Gamma_inv[t] = np.dot(deltaY[t - 1, :] / sigma2_e, weight.T) \
                + np.dot(mu_Gamma_inv[t - 1]
                         - deltaY_sigma2inv_rho_weightT[t - 1, :], tmp)
            mu[t] = np.linalg.solve(Gamma_inv[t], mu_Gamma_inv[t])
            tmp2 = mu_Gamma_inv[t - 1] - deltaY_sigma2inv_rho_weightT[t - 1, :]
            log_p_data += -self._half_log_det(Gamma_tilde_inv) \
                + np.dot(tmp2, np.linalg.solve(Gamma_tilde_inv, tmp2)) / 2.0
        log_p_data += -self._half_log_det(Gamma_inv[-1]) \
            + np.dot(mu_Gamma_inv[-1], mu[-1]) / 2.0 \
            - np.sum(deltaY**2 / sigma2_e) / 2.0
        return mu, mu_Gamma_inv, Gamma_inv, log_p_data, Lambda_0, \
            Lambda_1, H, deltaY, deltaY_sigma2inv_rho_weightT
1048578
def dstep(self, W):
        

        # Compute residual X D - S in frequency domain
        Ryf = sl.inner(self.Zf, self.Df, axis=self.cri.axisM) - self.Sf
        # Transform to spatial domain, apply mask, and transform back to
        # frequency domain
        Ryf[:] = sl.rfftn(W * sl.irfftn(Ryf, self.cri.Nv, self.cri.axisN),
                          None, self.cri.axisN)
        # Compute gradient
        gradf = sl.inner(np.conj(self.Zf), Ryf, axis=self.cri.axisK)

        # If multiple channel signal, single channel dictionary
        if self.cri.C > 1 and self.cri.Cd == 1:
            gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True)

        # Update gradient step
        self.eta = self.eta_a / (self.j + self.eta_b)

        # Compute gradient descent
        self.Gf[:] = self.Df - self.eta * gradf
        self.G = sl.irfftn(self.Gf, self.cri.Nv, self.cri.axisN)

        # Eval proximal operator
        self.Dprv[:] = self.D
        self.D[:] = self.Pcn(self.G)
999906
def conv1d(ni:int, no:int, ks:int=1, stride:int=1, padding:int=0, bias:bool=False):
    ""
    conv = nn.Conv1d(ni, no, ks, stride=stride, padding=padding, bias=bias)
    nn.init.kaiming_normal_(conv.weight)
    if bias: conv.bias.data.zero_()
    return spectral_norm(conv)
939687
def G(self):
        
        if self.c is None:
            raise RuntimeError('Inverse cross spectral density requires '
                               'invertible noise covariance matrix c.')
        A = self.A()
        # TODO: can we do that more efficiently?
        G = np.einsum('ji..., jk... ->ik...', A.conj(), self.Cinv())
        G = np.einsum('ij..., jk... ->ik...', G, A)
        return G
989980
def _dcorr(y, n2, A, dcov2_xx):
    
    # Pairwise Euclidean distances
    b = squareform(pdist(y, metric='euclidean'))
    # Double centering
    B = b - b.mean(axis=0)[None, :] - b.mean(axis=1)[:, None] + b.mean()
    # Compute squared distance covariances
    dcov2_yy = np.vdot(B, B) / n2
    dcov2_xy = np.vdot(A, B) / n2
    return np.sqrt(dcov2_xy) / np.sqrt(np.sqrt(dcov2_xx) * np.sqrt(dcov2_yy))
619129
def inference(self, kern, X, likelihood, Y, Y_metadata=None):

        
        N = X.shape[0] #number of training points
        D = X.shape[1] #number of dimensions

        Kds = np.zeros(D, dtype=object) #vector for holding covariance per dimension
        Qs = np.zeros(D, dtype=object) #vector for holding eigenvectors of covariance per dimension
        QTs = np.zeros(D, dtype=object) #vector for holding transposed eigenvectors of covariance per dimension
        V_kron = 1 # kronecker product of eigenvalues

        # retrieve the one-dimensional variation of the designated kernel
        oneDkernel = kern.get_one_dimensional_kernel(D)

        for d in range(D):
            xg = list(set(X[:,d])) #extract unique values for a dimension
            xg = np.reshape(xg, (len(xg), 1))
            oneDkernel.lengthscale = kern.lengthscale[d]
            Kds[d] = oneDkernel.K(xg)
            [V, Q] = np.linalg.eig(Kds[d])
            V_kron = np.kron(V_kron, V)
            Qs[d] = Q
            QTs[d] = Q.T

        noise = likelihood.variance + 1e-8

        alpha_kron = self.kron_mvprod(QTs, Y)
        V_kron = V_kron.reshape(-1, 1)
        alpha_kron = alpha_kron / (V_kron + noise)
        alpha_kron = self.kron_mvprod(Qs, alpha_kron)

        log_likelihood = -0.5 * (np.dot(Y.T, alpha_kron) + np.sum((np.log(V_kron + noise))) + N*log_2_pi)

        # compute derivatives wrt parameters Thete
        derivs = np.zeros(D+2, dtype='object')
        for t in range(len(derivs)):
            dKd_dTheta = np.zeros(D, dtype='object')
            gamma = np.zeros(D, dtype='object')
            gam = 1
            for d in range(D):
                xg = list(set(X[:,d]))
                xg = np.reshape(xg, (len(xg), 1))
                oneDkernel.lengthscale = kern.lengthscale[d]
                if t < D:
                    dKd_dTheta[d] = oneDkernel.dKd_dLen(xg, (t==d), lengthscale=kern.lengthscale[t]) #derivative wrt lengthscale
                elif (t == D):
                    dKd_dTheta[d] = oneDkernel.dKd_dVar(xg) #derivative wrt variance
                else:
                    dKd_dTheta[d] = np.identity(len(xg)) #derivative wrt noise
                gamma[d] = np.diag(np.dot(np.dot(QTs[d], dKd_dTheta[d].T), Qs[d]))
                gam = np.kron(gam, gamma[d])
            
            gam = gam.reshape(-1,1)
            kappa = self.kron_mvprod(dKd_dTheta, alpha_kron)
            derivs[t] = 0.5*np.dot(alpha_kron.T,kappa) - 0.5*np.sum(gam / (V_kron + noise))

        # separate derivatives
        dL_dLen = derivs[:D]
        dL_dVar = derivs[D]
        dL_dThetaL = derivs[D+1]

        return GridPosterior(alpha_kron=alpha_kron, QTs=QTs, Qs=Qs, V_kron=V_kron), \
                log_likelihood, {'dL_dLen':dL_dLen, 'dL_dVar':dL_dVar, 'dL_dthetaL':dL_dThetaL}
1269821
def S(self):
        
        if self.c is None:
            raise RuntimeError('Cross-spectral density requires noise '
                               'covariance matrix c.')
        H = self.H()
        # TODO: can we do that more efficiently?
        S = np.empty(H.shape, dtype=H.dtype)
        for f in range(H.shape[2]):
            S[:, :, f] = H[:, :, f].dot(self.c).dot(H[:, :, f].conj().T)
        return S
726356
def doubleprox_dc_simple(x, y, f, phi, g, K, niter, gamma, mu):
    
    for _ in range(niter):
        f.proximal(gamma)(x + gamma * K.adjoint(y) -
                          gamma * phi.gradient(x), out=x)
        g.convex_conj.proximal(mu)(y + mu * K(x), out=y)
823300
def gradient(self):
        """"""
        return (self.functional.gradient *
                (IdentityOperator(self.domain) - self.translation))