我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.diag()。
def annopred_inf(beta_hats, pr_sigi, n=1000, reference_ld_mats=None, ld_window_size=100): """ infinitesimal model with snp-specific heritability derived from annotation used as the initial values for MCMC of non-infinitesimal model """ num_betas = len(beta_hats) updated_betas = sp.empty(num_betas) m = len(beta_hats) for i, wi in enumerate(range(0, num_betas, ld_window_size)): start_i = wi stop_i = min(num_betas, wi + ld_window_size) curr_window_size = stop_i - start_i Li = 1.0/pr_sigi[start_i: stop_i] D = reference_ld_mats[i] A = (n/(1))*D + sp.diag(Li) A_inv = linalg.pinv(A) updated_betas[start_i: stop_i] = sp.dot(A_inv / (1.0/n) , beta_hats[start_i: stop_i]) # Adjust the beta_hats return updated_betas
def rhoA(self): # rhoA rhoA = pd.DataFrame(0, index=np.arange(1), columns=self.latent) for i in range(self.lenlatent): weights = pd.DataFrame(self.outer_weights[self.latent[i]]) weights = weights[(weights.T != 0).any()] result = pd.DataFrame.dot(weights.T, weights) result_ = pd.DataFrame.dot(weights, weights.T) S = self.data_[self.Variables['measurement'][ self.Variables['latent'] == self.latent[i]]] S = pd.DataFrame.dot(S.T, S) / S.shape[0] numerador = ( np.dot(np.dot(weights.T, (S - np.diag(np.diag(S)))), weights)) denominador = ( (np.dot(np.dot(weights.T, (result_ - np.diag(np.diag(result_)))), weights))) rhoA_ = ((result)**2) * (numerador / denominador) if(np.isnan(rhoA_.values)): rhoA[self.latent[i]] = 1 else: rhoA[self.latent[i]] = rhoA_.values return rhoA.T
def alpha(self): # Cronbach Alpha alpha = pd.DataFrame(0, index=np.arange(1), columns=self.latent) for i in range(self.lenlatent): block = self.data_[self.Variables['measurement'] [self.Variables['latent'] == self.latent[i]]] p = len(block.columns) if(p != 1): p_ = len(block) correction = np.sqrt((p_ - 1) / p_) soma = np.var(np.sum(block, axis=1)) cor_ = pd.DataFrame.corr(block) denominador = soma * correction**2 numerador = 2 * np.sum(np.tril(cor_) - np.diag(np.diag(cor_))) alpha_ = (numerador / denominador) * (p / (p - 1)) alpha[self.latent[i]] = alpha_ else: alpha[self.latent[i]] = 1 return alpha.T
def _compute_process_and_covariance_matrices(self, dt): """Computes the transition and covariance matrix of the process model and measurement model. Args: dt (float): Timestep of the discrete transition. Returns: F (numpy.ndarray): Transition matrix. Q (numpy.ndarray): Process covariance matrix. R (numpy.ndarray): Measurement covariance matrix. """ F = np.array(np.bmat([[np.eye(3), dt * np.eye(3)], [np.zeros((3, 3)), np.eye(3)]])) self.process_matrix = F q_p = self.process_covariance_position q_v = self.process_covariance_velocity Q = np.diag([q_p, q_p, q_p, q_v, q_v, q_v]) ** 2 * dt r = self.measurement_covariance R = r * np.eye(4) self.process_covariance = Q self.measurement_covariance = R return F, Q, R
def get_covariance(self): """Compute data covariance with the generative model. ``cov = components_.T * S**2 * components_ + sigma2 * eye(n_features)`` where S**2 contains the explained variances. Returns ------- cov : array, shape=(n_features, n_features) Estimated covariance of data. """ components_ = self.components_ exp_var = self.explained_variance_ if self.whiten: components_ = components_ * np.sqrt(exp_var[:, np.newaxis]) exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.) cov = np.dot(components_.T * exp_var_diff, components_) cov.flat[::len(cov) + 1] += self.noise_variance_ # modify diag inplace return cov
def get_scores(self): """Returns accuracy score evaluation result. - overall accuracy - mean accuracy - mean IU - fwavacc """ hist = self.confusion_matrix acc = np.diag(hist).sum() / hist.sum() acc_cls = np.diag(hist) / hist.sum(axis=1) acc_cls = np.nanmean(acc_cls) iu = np.diag(hist) / (hist.sum(axis=1) + hist.sum(axis=0) - np.diag(hist)) mean_iu = np.nanmean(iu) freq = hist.sum(axis=1) / hist.sum() fwavacc = (freq[freq > 0] * iu[freq > 0]).sum() cls_iu = dict(zip(range(self.n_classes), iu)) return {'Overall Acc: \t': acc, 'Mean Acc : \t': acc_cls, 'FreqW Acc : \t': fwavacc, 'Mean IoU : \t': mean_iu,}, cls_iu
def scale_variance(Theta, eps): """Allows to scale a Precision Matrix such that its corresponding covariance has unit variance Parameters ---------- Theta: ndarray Precision Matrix eps: float values to threshold to zero Returns ------- Theta: ndarray Precision of rescaled Sigma Sigma: ndarray Sigma with ones on diagonal """ Sigma = np.linalg.inv(Theta) V = np.diag(np.sqrt(np.diag(Sigma) ** -1)) Sigma = V.dot(Sigma).dot(V.T) # = VSV Theta = np.linalg.inv(Sigma) Theta[np.abs(Theta) <= eps] = 0. return Theta, Sigma
def sampson_error(F, pts1, pts2): """ Computes the sampson error for F, and points pts1, pts2. Sampson error is the first order approximation to the geometric error. Remember that this is a squared error. (x'^{T} * F * x)^2 ----------------- (F * x)_1^2 + (F * x)_2^2 + (F^T * x')_1^2 + (F^T * x')_2^2 where (F * x)_i^2 is the square of the i-th entry of the vector Fx """ x1, x2 = unproject_points(pts1).T, unproject_points(pts2).T Fx1 = np.dot(F, x1) Fx2 = np.dot(F, x2) # Sampson distance as error measure denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2 return ( np.diag(x1.T.dot(Fx2)) )**2 / denom
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = compute_xopt(self.rseed, dim) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = (self.condition ** .5) ** linspace(0, 1, dim) self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales)) # decouple scaling from function definition self.linearTF = dot(self.linearTF, self.rotation) # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape self.arrxopt = resize(self.xopt, curshape)
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = compute_xopt(self.rseed, dim) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = self.condition ** linspace(0, 1, dim) self.linearTF = dot(compute_rotation(self.rseed, dim), diag(((self.condition / 10.)**.5) ** linspace(0, 1, dim))) # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape self.arrxopt = resize(self.xopt, curshape)
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = compute_xopt(self.rseed, dim) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = (self.condition ** .5) ** linspace(0, 1, dim) self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales)) self.linearTF = dot(self.linearTF, self.rotation) # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape self.arrxopt = resize(self.xopt, curshape)
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = compute_xopt(self.rseed, dim) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = (self.condition ** .5) ** linspace(0, 1, dim) self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales)) # decouple scaling from function definition self.linearTF = dot(self.linearTF, self.rotation) # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape self.arrxopt = resize(self.xopt, curshape) self.arrexpo = resize(self.beta * linspace(0, 1, dim), curshape)
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = compute_xopt(self.rseed, dim) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = (1. / self.condition ** .5) ** linspace(0, 1, dim) # CAVE? self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales)) # decouple scaling from function definition self.linearTF = dot(self.linearTF, self.rotation) K = np.arange(0, 12) self.aK = np.reshape(0.5 ** K, (1, 12)) self.bK = np.reshape(3. ** K, (1, 12)) self.f0 = np.sum(self.aK * np.cos(2 * np.pi * self.bK * 0.5)) # optimal value # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape self.arrxopt = resize(self.xopt, curshape)
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = .5 * self._mu1 * sign(gauss(dim, self.rseed)) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = (self.condition ** .5) ** linspace(0, 1, dim) self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales)) # decouple scaling from function definition self.linearTF = dot(self.linearTF, self.rotation) # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape # self.arrxopt = resize(self.xopt, curshape) self.arrscales = resize(2. * sign(self.xopt), curshape) # makes up for xopt
def quaternion_matrix(quaternion): """Return homogeneous rotation matrix from quaternion. >>> M = quaternion_matrix([0.99810947, 0.06146124, 0, 0]) >>> numpy.allclose(M, rotation_matrix(0.123, [1, 0, 0])) True >>> M = quaternion_matrix([1, 0, 0, 0]) >>> numpy.allclose(M, numpy.identity(4)) True >>> M = quaternion_matrix([0, 1, 0, 0]) >>> numpy.allclose(M, numpy.diag([1, -1, -1, 1])) True """ q = numpy.array(quaternion, dtype=numpy.float64, copy=True) n = numpy.dot(q, q) if n < _EPS: return numpy.identity(4) q *= math.sqrt(2.0 / n) q = numpy.outer(q, q) return numpy.array([ [1.0-q[2, 2]-q[3, 3], q[1, 2]-q[3, 0], q[1, 3]+q[2, 0], 0.0], [ q[1, 2]+q[3, 0], 1.0-q[1, 1]-q[3, 3], q[2, 3]-q[1, 0], 0.0], [ q[1, 3]-q[2, 0], q[2, 3]+q[1, 0], 1.0-q[1, 1]-q[2, 2], 0.0], [ 0.0, 0.0, 0.0, 1.0]])
def distribution_parameters(self, parameter_name): if parameter_name=='trend': E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix) mean = dot(inv(eye(self.size)+E),self.data) cov = (self.parameters.list['sigma2'].current_value)*inv(eye(self.size)+E) return {'mean' : mean, 'cov' : cov} elif parameter_name=='sigma2': E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix) pos = self.size loc = 0 scale = 0.5*dot((self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)).T,(self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)))+0.5*dot(dot(self.parameters.list['trend'].current_value.T,E),self.parameters.list['trend'].current_value) elif parameter_name=='lambda2': pos = self.size-self.total_variation_order-1+self.alpha loc = 0.5*(norm(dot(self.derivative_matrix,self.parameters.list['trend'].current_value),ord=1))/self.parameters.list['sigma2'].current_value+self.rho scale = 1 elif parameter_name==str('omega'): pos = [sqrt(((self.parameters.list['lambda2'].current_value**2)*self.parameters.list['sigma2'].current_value)/(dj**2)) for dj in dot(self.derivative_matrix,self.parameters.list['trend'].current_value)] loc = 0 scale = self.parameters.list['lambda2'].current_value**2 return {'pos' : pos, 'loc' : loc, 'scale' : scale}
def test_psi(adjcube): """Tests retrieval of the wave functions and eigenvalues. """ from pydft.bases.fourier import psi, O, H cell = adjcube V = QHO(cell) W = W4(cell) Ns = W.shape[1] Psi, epsilon = psi(V, W, cell, forceR=False) #Make sure that the eigenvalues are real. assert np.sum(np.imag(epsilon)) < 1e-13 checkI = np.dot(Psi.conjugate().T, O(Psi, cell)) assert abs(np.sum(np.diag(checkI))-Ns) < 1e-13 # Should be the identity assert np.abs(np.sum(checkI)-Ns) < 1e-13 checkD = np.dot(Psi.conjugate().T, H(V, Psi, cell)) diagsum = np.sum(np.diag(checkD)) assert np.abs(np.sum(checkD)-diagsum) < 1e-12 # Should be diagonal # Should match the diagonal elements of previous matrix assert np.allclose(np.diag(checkD), epsilon)
def computePolarVecs(self,karg=False): N = len(self.times) L = np.reshape(self.L,(3,N)) if karg is False: A = self.computeRotMatrix() elif np.size(karg) is 3: A = self.computeRotMatrix(karg) elif np.size(karg) is 9: A = karg q = np.zeros((6,N)) for pp in range(0,N): Lpp = np.diag(L[:,pp]) p = np.dot(A,np.dot(Lpp,A.T)) q[:,pp] = np.r_[p[:,0],p[[1,2],1],p[2,2]] return q
def toa_normalize(x0, y0): xdim = x0.shape[0] m = x0.shape[1] n = x0.shape[1] t = -x0[:, 1] x = x0 + np.tile(t, (1, m)) y = y0 + np.tile(t, (1, n)) qr_a = x[:, 2:(1 + xdim)] q, r = scipy.linalg.qr(qr_a) x = (q.conj().T) * x y = (q.conj().T) * y M = np.diag(sign(np.diag(qr_a))) x1 = M * x y1 = M * y return x1, y1
def update_per_row(self, y_i, phi_i, J, mu0, c, v, r_prev_i, u_prev_i, phi_r_i, phi_u): # Params: # J - column indices nnz_i = len(J) residual_i = y_i - mu0 - c[J] prior_Phi = np.diag(np.concatenate(([phi_r_i], phi_u))) v_T = np.hstack((np.ones((nnz_i, 1)), v[J, :])) post_Phi_i = prior_Phi + \ np.dot(v_T.T, np.tile(phi_i[:, np.newaxis], (1, 1 + self.num_factor)) * v_T) # Weighted sum of v_j * v_j.T post_mean_i = np.squeeze(np.dot(phi_i * residual_i, v_T)) C, lower = scipy.linalg.cho_factor(post_Phi_i) post_mean_i = scipy.linalg.cho_solve((C, lower), post_mean_i) # Generate Gaussian, recycling the Cholesky factorization from the posterior mean computation. ru_i = math.sqrt(1 - self.relaxation ** 2) * scipy.linalg.solve_triangular(C, np.random.randn(len(post_mean_i)), lower=lower) ru_i += post_mean_i + self.relaxation * (post_mean_i - np.concatenate(([r_prev_i], u_prev_i))) r_i = ru_i[0] u_i = ru_i[1:] return r_i, u_i
def update_per_col(self, y_j, phi_j, I, mu0, r, u, c_prev_j, v_prev_j, phi_c_j, phi_v): prior_Phi = np.diag(np.concatenate(([phi_c_j], phi_v))) nnz_j = len(I) residual_j = y_j - mu0 - r[I] u_T = np.hstack((np.ones((nnz_j, 1)), u[I, :])) post_Phi_j = prior_Phi + \ np.dot(u_T.T, np.tile(phi_j[:, np.newaxis], (1, 1 + self.num_factor)) * u_T) # Weighted sum of u_i * u_i.T post_mean_j = np.squeeze(np.dot(phi_j * residual_j, u_T)) C, lower = scipy.linalg.cho_factor(post_Phi_j) post_mean_j = scipy.linalg.cho_solve((C, lower), post_mean_j) # Generate Gaussian, recycling the Cholesky factorization from the posterior mean computation. cv_j = math.sqrt(1 - self.relaxation ** 2) * scipy.linalg.solve_triangular(C, np.random.randn(len(post_mean_j)), lower=lower) cv_j += post_mean_j + self.relaxation * (post_mean_j - np.concatenate(([c_prev_j], v_prev_j))) c_j = cv_j[0] v_j = cv_j[1:] return c_j, v_j
def diag(v, k=0): """ Extract a diagonal or construct a diagonal array. This function is the equivalent of `numpy.diag` that takes masked values into account, see `numpy.diag` for details. See Also -------- numpy.diag : Equivalent function for ndarrays. """ output = np.diag(v, k).view(MaskedArray) if getmask(v) is not nomask: output._mask = np.diag(v._mask, k) return output
def predict_log_likelihood(self, paths, latents): if self.recurrent: observations = np.array([p["observations"][:, self.obs_regressed] for p in paths]) actions = np.array([p["actions"][:, self.act_regressed] for p in paths]) obs_actions = np.concatenate([observations, actions], axis=2) # latents must match first 2dim: (batch,time) else: observations = np.concatenate([p["observations"][:, self.obs_regressed] for p in paths]) actions = np.concatenate([p["actions"][:, self.act_regressed] for p in paths]) obs_actions = np.concatenate([observations, actions], axis=1) latents = np.concatenate(latents, axis=0) if self.noisify_traj_coef: noise = np.random.multivariate_normal(mean=np.zeros_like(np.mean(obs_actions, axis=0)), cov=np.diag(np.mean(np.abs(obs_actions), axis=0) * self.noisify_traj_coef), size=np.shape(obs_actions)[0]) obs_actions += noise if self.use_only_sign: obs_actions = np.sign(obs_actions) return self._regressor.predict_log_likelihood(obs_actions, latents) # see difference with fit above...
def lowb_mutual(self, paths, times=(0, None)): if self.recurrent: observations = np.array([p["observations"][times[0]:times[1], self.obs_regressed] for p in paths]) actions = np.array([p["actions"][times[0]:times[1], self.act_regressed] for p in paths]) obs_actions = np.concatenate([observations, actions], axis=2) latents = np.array([p['agent_infos']['latents'][times[0]:times[1]] for p in paths]) else: observations = np.concatenate([p["observations"][times[0]:times[1], self.obs_regressed] for p in paths]) actions = np.concatenate([p["actions"][times[0]:times[1], self.act_regressed] for p in paths]) obs_actions = np.concatenate([observations, actions], axis=1) latents = np.concatenate([p['agent_infos']["latents"][times[0]:times[1]] for p in paths]) if self.noisify_traj_coef: obs_actions += np.random.multivariate_normal(mean=np.zeros_like(np.mean(obs_actions,axis=0)), cov=np.diag(np.mean(np.abs(obs_actions), axis=0) * self.noisify_traj_coef), size=np.shape(obs_actions)[0]) if self.use_only_sign: obs_actions = np.sign(obs_actions) H_latent = self.policy.latent_dist.entropy(self.policy.latent_dist_info) # sum of entropies latents in return H_latent + np.mean(self._regressor.predict_log_likelihood(obs_actions, latents))
def gaussianWeight(kernelSize, even=False): if even == True: weight = np.ones([kernelSize,kernelSize]) weight = weight.reshape((1,kernelSize**2)) weight = np.array(weight)[0] weight = np.diag(weight) return weight SIGMA = 1 #the standard deviation of your normal curve CORRELATION = 0 #see wiki for multivariate normal distributions weight = np.zeros([kernelSize,kernelSize]) cpt = kernelSize%2 + kernelSize//2 #gets the center point for i in range(len(weight)): for j in range(len(weight)): ptx = i + 1 pty = j + 1 weight[i,j] = 1 / (2*np.pi*SIGMA**2) / (1-CORRELATION**2)**.5*np.exp(-1/(2*(1-CORRELATION**2))*((ptx-cpt)**2+(pty-cpt)**2)/(SIGMA**2)) # weight[i,j] = 1/SIGMA/(2*np.pi)**.5*np.exp(-(pt-cpt)**2/(2*SIGMA**2)) weight = weight.reshape((1,kernelSize**2)) weight = np.array(weight)[0] #convert to a 1D array weight = np.diag(weight) #convert to n**2xn**2 diagonal matrix return weight # return np.diag(weight)
def get_sigma(self): """Returns the co-variance matrix of the spline Return: (np.ndarray): Co-variance matrix for the EOS shape is (nxn) where n is the dof of the model """ sigma = self.get_option('spline_sigma') return np.diag((sigma * self.prior.get_dof())**2) # alpha = 3/self.shape()\ # * np.log(self.get_option('spline_max') # /self.get_option('spline_min')) # beta = np.log(1 + self.get_option('spline_sigma')) #return self.gram
def test_model_pq(self): """Test the model PQ matrix generation """ new_model = self.bayes.update(models={ 'simp': self.model.update_dof([4, 2])}) P, q = new_model._get_model_pq() epsilon = np.array([2, 1]) sigma = inv(np.diag(np.ones(2))) P_true = sigma q_true = -np.dot(epsilon, sigma) npt.assert_array_almost_equal(P, P_true, decimal=8) npt.assert_array_almost_equal(q, q_true, decimal=8)
def get_sigma(self, models): r"""Returns the variance matrix variance is .. math:: \Sigma_i = \sigma_t^2 + \frac{\sigma_x^2}{V_{CJ}} **Where** - :math:`\sigma_t` is the error in time measurements - :math:`\sigma_x` is the error in sensor position - :math:`V_{CJ}` is the detonation velocity see :py:meth:`F_UNCLE.Utils.Experiment.Experiment.get_sigma` """ eos = self.check_models(models)[0] return np.diag(np.ones(self.shape()))
def supcell(file, scale, outmode): from ababe.io.io import GeneralIO import os import numpy as np basefname = os.path.basename(file) gcell = GeneralIO.from_file(file) scale_matrix = np.diag(np.array(scale)) sc = gcell.supercell(scale_matrix) out = GeneralIO(sc) print("PROCESSING: {:}".format(file)) if outmode == 'stdio': out.write_file(fname=None, fmt='vasp') else: ofname = "{:}_SUPC.{:}".format(basefname.split('.')[0], outmode) out.write_file(ofname)
def do_seg_tests(net, iter, save_format, dataset, layer='score', gt='label'): n_cl = net.blobs[layer].channels if save_format: save_format = save_format.format(iter) hist, loss = compute_hist(net, save_format, dataset, layer, gt) # mean loss print '>>>', datetime.now(), 'Iteration', iter, 'loss', loss # overall accuracy acc = np.diag(hist).sum() / hist.sum() print '>>>', datetime.now(), 'Iteration', iter, 'overall accuracy', acc # per-class accuracy acc = np.diag(hist) / hist.sum(1) print '>>>', datetime.now(), 'Iteration', iter, 'mean accuracy', np.nanmean(acc) # per-class IU iu = np.diag(hist) / (hist.sum(1) + hist.sum(0) - np.diag(hist)) print '>>>', datetime.now(), 'Iteration', iter, 'mean IU', np.nanmean(iu) freq = hist.sum(1) / hist.sum() print '>>>', datetime.now(), 'Iteration', iter, 'fwavacc', \ (freq[freq > 0] * iu[freq > 0]).sum() return hist
def label_accuracy_score(label_trues, label_preds, n_class): """Returns accuracy score evaluation result. - overall accuracy - mean accuracy - mean IU - fwavacc """ hist = np.zeros((n_class, n_class)) for lt, lp in zip(label_trues, label_preds): hist += _fast_hist(lt.flatten(), lp.flatten(), n_class) acc = np.diag(hist).sum() / hist.sum() acc_cls = np.diag(hist) / hist.sum(axis=1) acc_cls = np.nanmean(acc_cls) iu = np.diag(hist) / (hist.sum(axis=1) + hist.sum(axis=0) - np.diag(hist)) mean_iu = np.nanmean(iu) freq = hist.sum(axis=1) / hist.sum() fwavacc = (freq[freq > 0] * iu[freq > 0]).sum() return acc, acc_cls, mean_iu, fwavacc # ----------------------------------------------------------------------------- # Visualization # -----------------------------------------------------------------------------
def fit(self, X): '''Required for featurewise_center, featurewise_std_normalization and zca_whitening. # Arguments X: Numpy array, the data to fit on. ''' if self.featurewise_center: self.mean = np.mean(X, axis=0) X -= self.mean if self.featurewise_std_normalization: self.std = np.std(X, axis=0) X /= (self.std + 1e-7) if self.zca_whitening: flatX = np.reshape(X, (X.shape[0], X.shape[1] * X.shape[2] * X.shape[3])) sigma = np.dot(flatX.T, flatX) / flatX.shape[1] U, S, V = linalg.svd(sigma) self.principal_components = np.dot(np.dot(U, np.diag(1. / np.sqrt(S + 10e-7))), U.T)
def __init__(self): """Initialize variable used by Kalman Filter class Args: None Return: None """ self.dt = 0.005 # delta time self.A = np.array([[1, 0], [0, 1]]) # matrix in observation equations self.u = np.zeros((2, 1)) # previous state vector # (x,y) tracking object center self.b = np.array([[0], [255]]) # vector of observations self.P = np.diag((3.0, 3.0)) # covariance matrix self.F = np.array([[1.0, self.dt], [0.0, 1.0]]) # state transition mat self.Q = np.eye(self.u.shape[0]) # process noise matrix self.R = np.eye(self.b.shape[0]) # observation noise matrix self.lastResult = np.array([[0], [255]])
def draw(vmean, vlogstd): from scipy import stats plt.cla() xlimits = [-2, 2] ylimits = [-4, 2] def log_prob(z): z1, z2 = z[:, 0], z[:, 1] return stats.norm.logpdf(z2, 0, 1.35) + \ stats.norm.logpdf(z1, 0, np.exp(z2)) plot_isocontours(ax, lambda z: np.exp(log_prob(z)), xlimits, ylimits) def variational_contour(z): return stats.multivariate_normal.pdf( z, vmean, np.diag(np.exp(vlogstd))) plot_isocontours(ax, variational_contour, xlimits, ylimits) plt.draw() plt.pause(1.0 / 30.0)
def get_neg_log_post(Phi, sigma_J_list, ROI_list, G, MMT, q, Sigma_E, GL, nu, V, prior_on = False): eps = 1E-13 p = Phi.shape[0] n_ROI = len(sigma_J_list) Qu = Phi.dot(Phi.T) G_Sigma_G = np.zeros(MMT.shape) for i in range(n_ROI): G_Sigma_G += sigma_J_list[i]**2 * np.dot(G[:,ROI_list[i]], G[:,ROI_list[i]].T) cov = Sigma_E + G_Sigma_G + GL.dot(Qu).dot(GL.T) inv_cov = np.linalg.inv(cov) eigs = np.real(np.linalg.eigvals(cov)) + eps log_det_cov = np.sum(np.log(eigs)) result = q*log_det_cov + np.trace(MMT.dot(inv_cov)) if prior_on: inv_Q = np.linalg.inv(Qu) #det_Q = np.linalg.det(Qu) log_det_Q = np.sum(np.log(np.diag(Phi)**2)) result = result + np.float(nu+p+1)*log_det_Q+ np.trace(V.dot(inv_Q)) return result #============================================================================== # update both Qu and Sigma_J, gradient of Qu and Sigma J
def fit(self, x): s = x.shape x = x.copy().reshape((s[0],np.prod(s[1:]))) m = np.mean(x, axis=0) x -= m sigma = np.dot(x.T,x) / x.shape[0] U, S, V = linalg.svd(sigma) tmp = np.dot(U, np.diag(1./np.sqrt(S+self.regularization))) tmp2 = np.dot(U, np.diag(np.sqrt(S+self.regularization))) self.ZCA_mat = th.shared(np.dot(tmp, U.T).astype(th.config.floatX)) self.inv_ZCA_mat = th.shared(np.dot(tmp2, U.T).astype(th.config.floatX)) self.mean = th.shared(m.astype(th.config.floatX))
def htmt(self): htmt_ = pd.DataFrame(pd.DataFrame.corr(self.data_), index=self.manifests, columns=self.manifests) mean = [] allBlocks = [] for i in range(self.lenlatent): block_ = self.Variables['measurement'][ self.Variables['latent'] == self.latent[i]] allBlocks.append(list(block_.values)) block = htmt_.ix[block_, block_] mean_ = (block - np.diag(np.diag(block))).values mean_[mean_ == 0] = np.nan mean.append(np.nanmean(mean_)) comb = [[k, j] for k in range(self.lenlatent) for j in range(self.lenlatent)] comb_ = [(np.sqrt(mean[comb[i][1]] * mean[comb[i][0]])) for i in range(self.lenlatent ** 2)] comb__ = [] for i in range(self.lenlatent ** 2): block = (htmt_.ix[allBlocks[comb[i][1]], allBlocks[comb[i][0]]]).values # block[block == 1] = np.nan comb__.append(np.nanmean(block)) htmt__ = np.divide(comb__, comb_) where_are_NaNs = np.isnan(htmt__) htmt__[where_are_NaNs] = 0 htmt = pd.DataFrame(np.tril(htmt__.reshape( (self.lenlatent, self.lenlatent)), k=-1), index=self.latent, columns=self.latent) return htmt
def _solve_equation_least_squares(self, A, B): """Solve system of linear equations A X = B. Currently using Pseudo-inverse because it also allows for singular matrices. Args: A (numpy.ndarray): Left-hand side of equation. B (numpy.ndarray): Right-hand side of equation. Returns: X (numpy.ndarray): Solution of equation. """ # Pseudo-inverse X = np.dot(np.linalg.pinv(A), B) # LU decomposition # lu, piv = scipy.linalg.lu_factor(A) # X = scipy.linalg.lu_solve((lu, piv), B) # Vanilla least-squares from numpy # X, _, _, _ = np.linalg.lstsq(A, B) # QR decomposition # Q, R, P = scipy.linalg.qr(A, mode='economic', pivoting=True) # # Find first zero element in R # out = np.where(np.diag(R) == 0)[0] # if out.size == 0: # i = R.shape[0] # else: # i = out[0] # B_prime = np.dot(Q.T, B) # X = np.zeros((A.shape[1], B.shape[1]), dtype=A.dtype) # X[P[:i], :] = scipy.linalg.solve_triangular(R[:i, :i], B_prime[:i, :]) return X
def get_whitening_matrix(X, fudge=1E-18): from numpy.linalg import eigh Xcov = numpy.dot(X.T, X)/X.shape[0] d,V = eigh(Xcov) D = numpy.diag(1./numpy.sqrt(d+fudge)) W = numpy.dot(numpy.dot(V,D), V.T) return W
def fit(self, X, C, y, regions, kernelType, reml=True, maxiter=100): #construct a list of kernel names (one for each region) if (kernelType == 'adapt'): kernelNames = self.buildKernelAdapt(X, C, y, regions, reml, maxiter) else: kernelNames = [kernelType] * len(regions) #perform optimization kernelObj, hyp_kernels, sig2e, fixedEffects = self.optimize(X, C, y, kernelNames, regions, reml, maxiter) #compute posterior distribution Ktraintrain = kernelObj.getTrainKernel(hyp_kernels) post = self.infExact_scipy_post(Ktraintrain, C, y, sig2e, fixedEffects) #fix intercept if phenotype is binary if (len(np.unique(y)) == 2): controls = (y<y.mean()) cases = ~controls meanVec = C.dot(fixedEffects) mu, var = self.getPosteriorMeanAndVar(np.diag(Ktraintrain), Ktraintrain, post, meanVec) fixedEffects[0] -= optimize.minimize_scalar(self.getNegLL, args=(mu, np.sqrt(sig2e+var), controls, cases), method='brent').x #construct trainObj trainObj = dict([]) trainObj['sig2e'] = sig2e trainObj['hyp_kernels'] = hyp_kernels trainObj['fixedEffects'] = fixedEffects trainObj['kernelNames'] = kernelNames return trainObj
def getPosteriorMeanAndVar(self, diagKTestTest, KtrainTest, post, intercept=0): L = post['L'] if (np.size(L) == 0): raise Exception('L is an empty array') #possible to compute it here Lchol = np.all((np.all(np.tril(L, -1)==0, axis=0) & (np.diag(L)>0)) & np.isreal(np.diag(L))) ns = diagKTestTest.shape[0] nperbatch = 5000 nact = 0 #allocate mem fmu = np.zeros(ns) #column vector (of length ns) of predictive latent means fs2 = np.zeros(ns) #column vector (of length ns) of predictive latent variances while (nact<(ns-1)): id = np.arange(nact, np.minimum(nact+nperbatch, ns)) kss = diagKTestTest[id] Ks = KtrainTest[:, id] if (len(post['alpha'].shape) == 1): try: Fmu = intercept[id] + Ks.T.dot(post['alpha']) except: Fmu = intercept + Ks.T.dot(post['alpha']) fmu[id] = Fmu else: try: Fmu = intercept[id][:, np.newaxis] + Ks.T.dot(post['alpha']) except: Fmu = intercept + Ks.T.dot(post['alpha']) fmu[id] = Fmu.mean(axis=1) if Lchol: V = la.solve_triangular(L, Ks*np.tile(post['sW'], (id.shape[0], 1)).T, trans=1, check_finite=False, overwrite_b=True) fs2[id] = kss - np.sum(V**2, axis=0) #predictive variances else: fs2[id] = kss + np.sum(Ks * (L.dot(Ks)), axis=0) #predictive variances fs2[id] = np.maximum(fs2[id],0) #remove numerical noise i.e. negative variances nact = id[-1] #set counter to index of last processed data point return fmu, fs2
def getTestKernelDiag(self, params, Xtest): self.checkParams(params) if (len(Xtest) != len(self.kernels)): raise Exception('Xtest should be a list with length equal to #kernels!') diag = 0 params_ind = 0 for k_i, k in enumerate(self.kernels): numHyp = k.getNumParams() diag += k.getTestKernelDiag(params[params_ind:params_ind+numHyp], Xtest[k_i]) params_ind += numHyp return diag #the only parameter here is the bias term c...
def getTestKernelDiag(self, params, Xtest): self.checkParams(params) b = np.exp(params[0]) k = np.exp(params[1]) M = 1.0 / (b + np.exp(k*self.D)) Xtest_scaled = Xtest/np.sqrt(Xtest.shape[1]) return np.diag((Xtest_scaled).dot(M).dot(Xtest_scaled.T)) #LD kernel with exponential decay but no bias term
def getTestKernelDiag(self, params, Xtest): self.checkParams(params) k = np.exp(params[0]) M = 1.0 / (1.0 + k*self.D) Xtest_scaled = Xtest/np.sqrt(Xtest.shape[1]) return np.diag((Xtest_scaled).dot(M).dot(Xtest_scaled.T)) #LD kernel with polynomial decay
def getTestKernelDiag(self, params, Xtest): self.checkParams(params) p = np.exp(params[0]) k = np.exp(params[1]) M = 1.0 / (1.0 + k*(self.D**p)) Xtest_scaled = Xtest/np.sqrt(Xtest.shape[1]) return np.diag((Xtest_scaled).dot(M).dot(Xtest_scaled.T))
def symmetrize(X): return X + X.T - np.diag(X.diagonal()) #this code is directly translated from the GPML Matlab package (see attached license)
def l1_norm_off_diag(A): "convenience method for l1 norm, excluding diagonal" # let's speed this up later # assume A symmetric, sparse too return np.linalg.norm(A - np.diag(A.diagonal()), ord=1)
