我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.cumsum()。
def reset_index(self): """Reset index to range based """ dfs = self.to_delayed() sizes = np.asarray(compute(*map(delayed(len), dfs))) prefixes = np.zeros_like(sizes) prefixes[1:] = np.cumsum(sizes[:-1]) @delayed def fix_index(df, startpos): return df.set_index(np.arange(start=startpos, stop=startpos + len(df), dtype=np.intp)) outdfs = [fix_index(df, startpos) for df, startpos in zip(dfs, prefixes)] return from_delayed(outdfs)
def scalar_weight_updates(xc, ec, kp, kd): tx_last = 0 te_last = 0 x_last = 0 e_last = 0 n_samples = xc.shape[0] dw = np.zeros(n_samples) r = kd/float(kp+kd) for t in xrange(n_samples): t_last = max(tx_last, te_last) if xc[t]: dw[t] = x_last*e_last*r**(2*t_last-tx_last-te_last)*geosum(r**2, t_end=t-t_last+1, t_start=1) x_last = x_last * r**(t-tx_last) + xc[t] tx_last = t if ec[t]: dw[t] = x_last*e_last*r**(2*t_last-tx_last-te_last)*geosum(r**2, t_end=t-t_last+1, t_start=1) # THing... e_last = e_last * r**(t-te_last) + ec[t] te_last = t return np.cumsum(dw)/kd**2
def pd_weight_grads_future(xc, ec, kp, kd): r = kd/float(kp+kd) kd=float(kd) scale_factor = 1./float(kp**2 + 2*kp*kd) n_samples = xc.shape[0] assert n_samples == ec.shape[0] n_in = xc.shape[1] n_out = ec.shape[1] dws = np.zeros((n_samples, n_in, n_out)) xr = np.zeros(n_in) er = np.zeros(n_out) for t in xrange(n_samples): xr *= r er = er*r + ec[t] dws[t] = (xc[t][:, None]*er[None, :] + xr[:, None]*ec[t][None, :]) * scale_factor xr += xc[t] # er += ec[t] return np.cumsum(dws, axis=0)
def _train_pca(self, training_set): """Trains and returns a LinearMachine that is trained using PCA""" data = numpy.vstack(feature for feature in training_set) logger.info(" -> Training LinearMachine using PCA ") trainer = bob.learn.linear.PCATrainer() machine, eigen_values = trainer.train(data) if isinstance(self.subspace_dimension_pca, float): cummulated = numpy.cumsum(eigen_values) / numpy.sum(eigen_values) for index in range(len(cummulated)): if cummulated[index] > self.subspace_dimension_pca: self.subspace_dimension_pca = index break self.subspace_dimension_pca = index # limit number of pcs logger.info(" -> limiting PCA subspace to %d dimensions", self.subspace_dimension_pca) machine.resize(machine.shape[0], self.subspace_dimension_pca) return machine
def train_projector(self, training_features, projector_file): """Generates the PCA covariance matrix and writes it into the given projector_file. **Parameters:** training_features : [1D :py:class:`numpy.ndarray`] A list of 1D training arrays (vectors) to train the PCA projection matrix with. projector_file : str A writable file, into which the PCA projection matrix (as a :py:class:`bob.learn.linear.Machine`) and the eigenvalues will be written. """ # Assure that all data are 1D [self._check_feature(feature) for feature in training_features] # Initializes the data data = numpy.vstack(training_features) logger.info(" -> Training LinearMachine using PCA") t = bob.learn.linear.PCATrainer() self.machine, self.variances = t.train(data) # For re-shaping, we need to copy... self.variances = self.variances.copy() # compute variance percentage, if desired if isinstance(self.subspace_dim, float): cummulated = numpy.cumsum(self.variances) / numpy.sum(self.variances) for index in range(len(cummulated)): if cummulated[index] > self.subspace_dim: self.subspace_dim = index break self.subspace_dim = index logger.info(" ... Keeping %d PCA dimensions", self.subspace_dim) # re-shape machine self.machine.resize(self.machine.shape[0], self.subspace_dim) self.variances.resize(self.subspace_dim) f = bob.io.base.HDF5File(projector_file, "w") f.set("Eigenvalues", self.variances) f.create_group("Machine") f.cd("/Machine") self.machine.save(f)
def split(x, split_dim, split_sizes): n = len(list(x.get_shape())) dim_size = np.sum(split_sizes) assert int(x.get_shape()[split_dim]) == dim_size ids = np.cumsum([0] + split_sizes) ids[-1] = -1 begin_ids = ids[:-1] ret = [] for i in range(len(split_sizes)): cur_begin = np.zeros([n], dtype=np.int32) cur_begin[split_dim] = begin_ids[i] cur_end = np.zeros([n], dtype=np.int32) - 1 cur_end[split_dim] = split_sizes[i] ret += [tf.slice(x, cur_begin, cur_end)] return ret
def GenCR(MCMCPar,pCR): if type(pCR) is np.ndarray: p=np.ndarray.tolist(pCR)[0] else: p=pCR CR=np.zeros((MCMCPar.seq * MCMCPar.steps),dtype=np.float) L = np.random.multinomial(MCMCPar.seq * MCMCPar.steps, p, size=1) L2 = np.concatenate((np.zeros((1),dtype=np.int), np.cumsum(L)),axis=0) r = np.random.permutation(MCMCPar.seq * MCMCPar.steps) for zz in range(0,MCMCPar.nCR): i_start = L2[zz] i_end = L2[zz+1] idx = r[i_start:i_end] CR[idx] = np.float(zz+1)/MCMCPar.nCR CR = np.reshape(CR,(MCMCPar.seq,MCMCPar.steps)) return CR, L
def DEStrategy(DEpairs,seq): # Determine which sequences to evolve with what DE strategy # Determine probability of selecting a given number of pairs p_pair = (1.0/DEpairs) * np.ones((1,DEpairs)) p_pair = np.cumsum(p_pair) p_pair = np.concatenate((np.zeros((1)),p_pair),axis=0) DEversion=np.zeros((seq),dtype=np.int32) Z = np.random.rand(seq) # Select number of pairs for qq in range(0,seq): z = np.where(p_pair<=Z[qq]) DEversion[qq] = z[0][-1] return DEversion
def cumultativesumstest(binin): ''' The focus of this test is the maximal excursion (from zero) of the random walk defined by the cumulative sum of adjusted (-1, +1) digits in the sequence. The purpose of the test is to determine whether the cumulative sum of the partial sequences occurring in the tested sequence is too large or too small relative to the expected behavior of that cumulative sum for random sequences. This cumulative sum may be considered as a random walk. For a random sequence, the random walk should be near zero. For non-random sequences, the excursions of this random walk away from zero will be too large.''' n = len(binin) ss = [int(el) for el in binin] sc = map(sumi, ss) cs = np.cumsum(sc) z = max(abs(cs)) ra = 0 start = int(np.floor(0.25 * np.floor(-n / z) + 1)) stop = int(np.floor(0.25 * np.floor(n / z) - 1)) pv1 = [] for k in xrange(start, stop + 1): pv1.append(sst.norm.cdf((4 * k + 1) * z / np.sqrt(n)) - sst.norm.cdf((4 * k - 1) * z / np.sqrt(n))) start = int(np.floor(0.25 * np.floor(-n / z - 3))) stop = int(np.floor(0.25 * np.floor(n / z) - 1)) pv2 = [] for k in xrange(start, stop + 1): pv2.append(sst.norm.cdf((4 * k + 3) * z / np.sqrt(n)) - sst.norm.cdf((4 * k + 1) * z / np.sqrt(n))) pval = 1 pval -= reduce(su, pv1) pval += reduce(su, pv2) return pval
def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" exType = exType.lower() possibleTypes = ['rect', 'rotate'] assert exType in possibleTypes, "Not a possible example type." if exType == 'rect': return list(ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)) elif exType == 'rotate': if len(nC) == 2: X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2) amt[amt < 0] = 0 return [X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt] elif len(nC) == 3: X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2) amt[amt < 0] = 0 return [X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt]
def test_basic(self): ba = [1, 2, 10, 11, 6, 5, 4] ba2 = [[1, 2, 3, 4], [5, 6, 7, 9], [10, 3, 4, 5]] for ctype in [np.int8, np.uint8, np.int16, np.uint16, np.int32, np.uint32, np.float32, np.float64, np.complex64, np.complex128]: a = np.array(ba, ctype) a2 = np.array(ba2, ctype) tgt = np.array([1, 3, 13, 24, 30, 35, 39], ctype) assert_array_equal(np.cumsum(a, axis=0), tgt) tgt = np.array( [[1, 2, 3, 4], [6, 8, 10, 13], [16, 11, 14, 18]], ctype) assert_array_equal(np.cumsum(a2, axis=0), tgt) tgt = np.array( [[1, 3, 6, 10], [5, 11, 18, 27], [10, 13, 17, 22]], ctype) assert_array_equal(np.cumsum(a2, axis=1), tgt)
def __init__(self, delimiter=None, comments=asbytes('#'), autostrip=True): self.comments = comments # Delimiter is a character if isinstance(delimiter, unicode): delimiter = delimiter.encode('ascii') if (delimiter is None) or _is_bytes_like(delimiter): delimiter = delimiter or None _handyman = self._delimited_splitter # Delimiter is a list of field widths elif hasattr(delimiter, '__iter__'): _handyman = self._variablewidth_splitter idx = np.cumsum([0] + list(delimiter)) delimiter = [slice(i, j) for (i, j) in zip(idx[:-1], idx[1:])] # Delimiter is a single integer elif int(delimiter): (_handyman, delimiter) = ( self._fixedwidth_splitter, int(delimiter)) else: (_handyman, delimiter) = (self._delimited_splitter, None) self.delimiter = delimiter if autostrip: self._handyman = self.autostrip(_handyman) else: self._handyman = _handyman #
def map_index(self, index): # Get a list of lengths of all datasets. Say the answer is [4, 3, 3], # and we're looking for index = 5. len_list = list(map(len, self.datasets)) # Cumulate to a numpy array. The answer is [4, 7, 10] cumulative_len_list = np.cumsum(len_list) # When the index is subtracted, we get [-1, 2, 5]. We're looking for the (index # of the) first cumulated len which is larger than the index (in this case, # 7 (index 1)). offset_cumulative_len_list = cumulative_len_list - index dataset_index = np.argmax(offset_cumulative_len_list > 0) # With the dataset index, we figure out the index in dataset if dataset_index == 0: # First dataset - index corresponds to index_in_dataset index_in_dataset = index else: # Get cumulated length up to the current dataset len_up_to_dataset = cumulative_len_list[dataset_index - 1] # Compute index_in_dataset as that what's left index_in_dataset = index - len_up_to_dataset return dataset_index, index_in_dataset
def cumsum(a, axis=None, dtype=None, out=None): """Returns the cumulative sum of an array along a given axis. Args: a (cupy.ndarray): Input array. axis (int): Axis along which the cumulative sum is taken. If it is not specified, the input is flattened. dtype: Data type specifier. out (cupy.ndarray): Output array. Returns: cupy.ndarray: The result array. .. seealso:: :func:`numpy.cumsum` """ return _cum_core(a, axis, dtype, out, _cumsum_kern, _cumsum_batch_kern)
def __call__(self, time_sequence, weather_data): """ Compute thermal time accumulation over time_sequence :Parameters: ---------- - `time_sequence` (panda dateTime index) A sequence of TimeStamps indicating the dates of all elementary time steps of the simulation - weather (alinea.astk.Weather instance) A Weather database """ try: Tair = weather_data.temperature_air[time_sequence] except: #strange extract needed on visualea 1.0 (to test again with ipython in visualea) T_data = weather_data[['temperature_air']] Tair = numpy.array([float(T_data.loc[d]) for d in time_sequence]) Tcut = numpy.maximum(numpy.zeros_like(Tair), Tair - self.Tbase) days = [0] + [((t - time_sequence[0]).total_seconds()+ 3600) / 3600 / 24 for t in time_sequence] dt = numpy.diff(days).tolist() return numpy.cumsum(Tcut * dt) # functional call for nodes
def accumulate_grad(self): if self.n_backward == 0: self.model.zerograds() # Compute losses losses = [] for r_seq, log_prob_seq, ent_seq in zip(self.reward_sequences, self.log_prob_sequences, self.entropy_sequences): assert len(r_seq) - 1 == len(log_prob_seq) == len(ent_seq) # Convert rewards into returns (=sum of future rewards) R_seq = np.cumsum(list(reversed(r_seq[1:])))[::-1] for R, log_prob, entropy in zip(R_seq, log_prob_seq, ent_seq): loss = -R * log_prob - self.beta * entropy losses.append(loss) total_loss = chainerrl.functions.sum_arrays(losses) # When self.batchsize is future.types.newint.newint, dividing a # Variable with it will raise an error, so it is manually converted to # float here. total_loss /= float(self.batchsize) total_loss.backward() self.reward_sequences = [[]] self.log_prob_sequences = [[]] self.entropy_sequences = [[]] self.n_backward += 1
def _to_notesequences(self, samples): output_sequences = [] dim_ranges = np.cumsum(self._split_output_depths) for s in samples: mel_ns = self._melody_converter.to_notesequences( [s[:, :dim_ranges[0]]])[0] bass_ns = self._melody_converter.to_notesequences( [s[:, dim_ranges[0]:dim_ranges[1]]])[0] drums_ns = self._drums_converter.to_notesequences( [s[:, dim_ranges[1]:]])[0] for n in bass_ns.notes: n.instrument = 1 n.program = ELECTRIC_BASS_PROGRAM for n in drums_ns.notes: n.instrument = 9 ns = mel_ns ns.notes.extend(bass_ns.notes) ns.notes.extend(drums_ns.notes) ns.total_time = max( mel_ns.total_time, bass_ns.total_time, drums_ns.total_time) output_sequences.append(ns) return output_sequences
def sample_categorical(pmf): """Sample from a categorical distribution. Args: pmf: Probablity mass function. Output of a softmax over categories. Array of shape [batch_size, number of categories]. Rows sum to 1. Returns: idxs: Array of size [batch_size, 1]. Integer of category sampled. """ if pmf.ndim == 1: pmf = np.expand_dims(pmf, 0) batch_size = pmf.shape[0] cdf = np.cumsum(pmf, axis=1) rand_vals = np.random.rand(batch_size) idxs = np.zeros([batch_size, 1]) for i in range(batch_size): idxs[i] = cdf[i].searchsorted(rand_vals[i]) return idxs
def simBirth(self,which_agents): ''' Makes new Markov consumer by drawing initial normalized assets, permanent income levels, and discrete states. Calls IndShockConsumerType.simBirth, then draws from initial Markov distribution. Parameters ---------- which_agents : np.array(Bool) Boolean array of size self.AgentCount indicating which agents should be "born". Returns ------- None ''' IndShockConsumerType.simBirth(self,which_agents) # Get initial assets and permanent income N = np.sum(which_agents) base_draws = drawUniform(N,seed=self.RNG.randint(0,2**31-1)) Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit)) self.MrkvNow[which_agents] = np.searchsorted(Cutoffs,base_draws).astype(int)
def clip_spectrum(s, k, discard=0.001): """ Given eigenvalues `s`, return how many factors should be kept to avoid storing spurious (tiny, numerically instable) values. This will ignore the tail of the spectrum with relative combined mass < min(`discard`, 1/k). The returned value is clipped against `k` (= never return more than `k`). """ # compute relative contribution of eigenvalues towards the energy spectrum rel_spectrum = numpy.abs(1.0 - numpy.cumsum(s / numpy.sum(s))) # ignore the last `discard` mass (or 1/k, whichever is smaller) of the spectrum small = 1 + len(numpy.where(rel_spectrum > min(discard, 1.0 / k))[0]) k = min(k, small) # clip against k logger.info("keeping %i factors (discarding %.3f%% of energy spectrum)" % (k, 100 * rel_spectrum[k - 1])) return k
def plot_stack_cum(semantic, *reports, color_theme = C_cold_colors, color_offset = 0): """ Creates a stacked plot with cumulated sums for a given semantic """ X, Y, c = extract_data(semantic, *reports, color_theme = color_theme, color_offset = color_offset) Y = remove_nones(Y) # add zeros only at the beginning for x in X: x.insert(0, x[0] - 1) for data in Y: for d in data: d.insert(0, 0.) # create the cumulated sum Y = [np.cumsum(np.array(yi), 1) if yi else [] for yi in Y] plot_stack_generic(X, Y, c, semantic, *reports, color_theme = color_theme, color_offset = color_offset) legend(loc = 'upper left', fancybox=True, framealpha=0.4, prop={'size':10})
def get_segment_vote_accuracy(segment_label, segment_predictions, window): def gen(): count = { label: np.hstack([[0], np.cumsum(segment_predictions == label)]) for label in set(segment_predictions) } tmp = window if tmp == -1: tmp = len(segment_predictions) tmp = min(tmp, len(segment_predictions)) for begin in range(len(segment_predictions) - tmp + 1): yield segment_label == max( count, key=lambda label: count[label][begin + tmp] - count[label][begin] ), segment_label return list(gen())
def _parse_headers(self): self.num_data_list = [] self.ones_accum_list = [] self.multi_accum_list = [] self.num_pix = [] for i, photons_file in enumerate(self.photons_list): with open(photons_file, 'rb') as f: num_data = np.fromfile(f, dtype='i4', count=1)[0] self.num_pix.append(np.fromfile(f, dtype='i4', count=1)[0]) if self.num_pix[i] != len(self.geom_list[i].x): sys.stderr.write('Warning: num_pix for %s is different (%d vs %d)\n' % (photons_file, self.num_pix[i], len(self.geom_list[i].x))) f.seek(1024, 0) ones = np.fromfile(f, dtype='i4', count=num_data) multi = np.fromfile(f, dtype='i4', count=num_data) self.num_data_list.append(num_data) self.ones_accum_list.append(np.cumsum(ones)) self.multi_accum_list.append(np.cumsum(multi)) self.num_data_list = np.cumsum(self.num_data_list) self.num_frames = self.num_data_list[-1]
def fit_koff(nmax=523, NN=4e8, **params): tbind = params.pop("tbind") params["kd"] = 1e9/tbind dx = params.pop("dx") rw = randomwalk.get_rw(NAME, params, setup=setup_rw, calc=True) rw.domains[1].dx = dx times = draw_empirically(rw, N=NN, nmax=nmax, success=False) bins = np.logspace(np.log10(min(times)), np.log10(max(times)), 35) #bins = np.logspace(-3., 2., 35) hist, _ = np.histogram(times, bins=bins) cfd = np.cumsum(hist)/float(np.sum(hist)) t = 0.5*(bins[:-1] + bins[1:]) tmean = times.mean() toff = NLS(t, cfd, t0=tmean) koff = 1./toff return dict(t=t, cfd=cfd, toff=toff, tmean=tmean, koff=koff) ##### run rw in collect mode and draw bindings from empirical distributions
def convert(batch, device): def to_device_batch(batch): if device is None: return batch elif device < 0: return [chainer.dataset.to_device(device, x) for x in batch] else: xp = cuda.cupy.get_array_module(*batch) concat = xp.concatenate(batch, axis=0) sections = np.cumsum([len(x) for x in batch[:-1]], dtype='i') concat_dev = chainer.dataset.to_device(device, concat) batch_dev = cuda.cupy.split(concat_dev, sections) return batch_dev return {'xs': to_device_batch([x for x, _ in batch]), 'ys': to_device_batch([y for _, y in batch])}
def csc_matvec(mat_csc, vec, dense_output=True, dtype=None): v_nnz = vec.indices v_val = vec.data m_val = mat_csc.data m_ind = mat_csc.indices m_ptr = mat_csc.indptr res_dtype = dtype or np.result_type(mat_csc.dtype, vec.dtype) if dense_output: res = np.zeros((mat_csc.shape[0],), dtype=res_dtype) matvec2dense(m_ptr, m_ind, m_val, v_nnz, v_val, res) else: sizes = m_ptr.take(v_nnz+1) - m_ptr.take(v_nnz) sizes = np.concatenate(([0], np.cumsum(sizes))) n = sizes[-1] data = np.empty((n,), dtype=res_dtype) indices = np.empty((n,), dtype=np.intp) indptr = np.array([0, n], dtype=np.intp) matvec2sparse(m_ptr, m_ind, m_val, v_nnz, v_val, sizes, indices, data) res = sp.sparse.csr_matrix((data, indices, indptr), shape=(1, mat_csc.shape[0]), dtype=res_dtype) res.sum_duplicates() # expensive operation return res
def plot_loss(loss_list, log_dir, iter_id): def running_mean(x, N): cumsum = np.cumsum(np.insert(x, 0, 0)) return (cumsum[N:] - cumsum[:-N]) / N plt.figure() plt.semilogy(loss_list, '.', alpha=0.2, label="Loss") plt.semilogy(running_mean(loss_list,100), label="Average Loss") plt.xlabel('Iterations') plt.ylabel('Loss') plt.legend() plt.grid() ax = plt.subplot(111) ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.05), ncol=3, fancybox=True, shadow=True) plt.savefig(log_dir + "/fig_loss_iter_" + str(iter_id) + ".pdf") print("figure plotted") plt.close()
def _moving_average(self, a): """Compute the moving average of a signal. Parameters ---------- a : np.array Signal Returns ------- a_avg : np.array Moving average of signal. """ n = self.avg_filt_len ret = np.cumsum(a, dtype=float) ret[n:] = ret[n:] - ret[:-n] return ret[n - 1:] / n
def PCAdo(block, name): cor_ = np.corrcoef(block.T) eig_vals, eig_vecs = np.linalg.eig(cor_) tot = sum(eig_vals) var_exp = [(i / tot) * 100 for i in sorted(eig_vals, reverse=True)] cum_var_exp = np.cumsum(var_exp) loadings = (eig_vecs * np.sqrt(eig_vals)) eig_vals = np.sort(eig_vals)[::-1] print('Eigenvalues') print(eig_vals) print('Variance Explained') print(var_exp) print('Total Variance Explained') print(cum_var_exp) print('Loadings') print(abs(loadings[:, 0])) PAcorrect = PA(block.shape[0], block.shape[1]) print('Parallel Analisys') pa = (eig_vals - (PAcorrect - 1)) print(pa) print('Correlation Matrix') print(pd.DataFrame.corr(block)) plt.plot(range(1,len(pa)+1), pa, '-o') plt.grid(True) plt.xlabel('Fatores') plt.ylabel('Componentes') plt.savefig('imgs/PCA' + name, bbox_inches='tight') plt.clf() plt.cla() # plt.show()
def init_ecdf(self): if self.count == 0: raise Exception('Need to add items before you can call init_ecdf()') self.init_prob() self._ecdf = np.cumsum( self._prob ) return # Precondition: must call init_prob() after last add() operation
def reconstruct(self): # reconstruct progressively self.reconstruction[:, self.output_ranks] = np.cumsum( self.weights[:, self.output_ranks] * self.output[:, self.output_ranks], axis=1) self.activation(self.reconstruction, out=self.reconstruction)
def get_sparsity_error_term(self): # the total error function being minimized return np.sum(np.cumsum(self.get_reconstruction_error_vector()[::-1]))
def moving_average(a, n=3): ret = np.cumsum(a) ret[n:] = ret[n:] - ret[:-n] return ret[n - 1:] / n
def fit_optimal(blocks, keys): """Fit the cumulative distribution, by optimization.""" values = list() for k in keys: values.extend([max_val for (_, _, max_val) in blocks[k]]) hist, bin_edges = np.histogram(values, bins=100, normed=True) bin_centers = np.array( [(bin_edges[i - 1] + bin_edges[i]) / 2 for i in range(1, len(bin_edges))]) cumulative = np.cumsum(hist) / np.sum(hist) popt, _ = optimize.curve_fit(func, bin_centers, cumulative) return popt
def pd_weight_grads(xc, ec, kp, kd): """ Efficiently compute the weights over time for a system recieving sparse inputs xc and ec. :param xc: An (n_samples, n_in) array of input spikes :param ec: An (n_samples, n_in) array of error_spikes :param kp: A scalar kp value :param kd: A scalar kd value :return: An (n_samples, n_in, n_out) array of weight updates """ r = kd/float(kp+kd) n_samples = xc.shape[0] assert n_samples == ec.shape[0] n_in = xc.shape[1] n_out = ec.shape[1] dws = np.zeros((n_samples, n_in, n_out)) tx_last = np.zeros(n_in) te_last = np.zeros(n_out) x_last = np.zeros(n_in) e_last = np.zeros(n_out) for t in xrange(n_samples): x_spikes = xc[t] != 0 t_last = np.maximum(tx_last[x_spikes, None], te_last) dws[t, x_spikes, :] = x_last[x_spikes, None] * e_last * r**(2*t_last-tx_last[x_spikes, None]-te_last)*geosum(r**2, t_end=t-t_last, t_start=1) # Not 100% on this... x_last[x_spikes] = x_last[x_spikes]*r**(t-tx_last[x_spikes]) + xc[t][x_spikes] / float(kd) tx_last[x_spikes] = t e_spikes = ec[t] != 0 if np.any(e_spikes): t_last = np.maximum(tx_last[:, None], te_last[e_spikes]) dws[t, :, e_spikes] += (x_last[:, None] * e_last[e_spikes] * r**(2*t_last-tx_last[:, None]-te_last[e_spikes])*geosum(r**2, t_end=t-t_last, t_start=1)).T # T makes no sense here but e_last[e_spikes] = e_last[e_spikes]*r**(t-te_last[e_spikes]) + ec[t][e_spikes] / float(kd) te_last[e_spikes] = t return np.cumsum(dws, axis=0)
def pd_weight_grads_mult(xc, ec, kp, kd): """ Efficiently compute the weights over time for a system recieving sparse inputs xc and ec. We use a slightly different approach this time - instead of keeping the last values and spike times for each neuron, we :param xc: An (n_samples, n_in) array of input spikes :param ec: An (n_samples, n_in) array of error_spikes :param kp: A scalar kp value :param kd: A scalar kd value :return: An (n_samples, n_in, n_out) array of weight updates """ # TODO: Make this work and pass The Test r = kd/float(kp+kd) kd=float(kd) n_samples = xc.shape[0] assert n_samples == ec.shape[0] n_in = xc.shape[1] n_out = ec.shape[1] dws = np.zeros((n_samples, n_in, n_out)) xr = np.zeros(n_in) xi = np.zeros(n_in) er = np.zeros(n_out) ei = np.zeros(n_out) for t in xrange(n_samples): xr = r*xr + xc[t]/kd**2 er = r*er + ec[t]/kd**2 xi = r*xi*(1-(ec[t]!=0)) + xr ei = r*ei*((1-xc[t]!=0)) + er # dws[t] = xc[t, :, None]*ei[None, :] + xi[:, None]*ec[t, None, :] dws[t] = xc[t, :, None]*ei[None, :] + xi[:, None]*ec[t, None, :] return np.cumsum(dws, axis=0)*kd
def sequential_quantize(v, n_steps = None, method='herd', rng = None): """ :param v: A (..., n_samples, n_units, ) array :param n_steps: The number of steps to spike for :return: An (..., n_steps, n_units) array of quantized values """ rng = get_rng(rng) assert v.ndim>=2 if n_steps is None: n_steps = v.shape[-2] else: assert n_steps == v.shape[-2] if method=='herd': result = fixed_diff(np.round(np.cumsum(v, axis=-2)), axis=-2) elif method=='herd2': result = fixed_diff(fixed_diff(np.round(np.cumsum(np.cumsum(v, axis=-2), axis=-2)), axis=-2), axis=-2) elif method=='round': result = np.round(v) elif method == 'slippery.9': result = slippery_round(v, slip=0.9) elif method == 'slippery.5': result = slippery_round(v, slip=0.5) elif method == 'randn': result = v + rng.randn(*v.shape) elif method=='uniform': result = v + rng.uniform(-.5, .5, size=v.shape) elif method=='surrogate-noise': result = v + (12**.5)*((v%1)-(v%1)**2)*rng.uniform(low=-.5, high=.5, size=v.shape) elif method == 'surrogate-sqrt': result = v + np.sqrt((12**.5)*((v%1)-(v%1)**2)*rng.uniform(low=-.5, high=.5, size=v.shape)) elif method is None: result = v else: raise NotImplementedError("Don't have quantization method '%s' implemented" % (method, )) return result
def pitch_strength_all_candidates(f_erbs, L, pc): """ Calculates the pitch ``strength'' of all candidate pitches Args: f_erbs (array): frequencies in ERBs L (matrix): loudness matrix pc (array): pitch candidates array Returns: S (array): strength of pitches corresponding to pc's """ # create pitch strength matrix S = np.zeros((pc.size, L.shape[1])) # define integration regions k = np.zeros(pc.size+1) for j in range(k.size-1): idx = int(k[j]) f = f_erbs[idx:] val = find(f > pc[j] / 4)[0] k[j+1] = k[j] + val k = k[1:] # TODO: fix this sloppiness # create loudness normalization matrix N = np.sqrt(np.flipud(np.cumsum(np.flipud(L * L), 0))) for j in range(pc.size): # normalize loudness n = N[int(k[j]), :] n[n == 0] = -np.inf # to make zero-loudness equal zero after normalization nL = L[int(k[j]):] / np.tile(n, (int(L.shape[0] - k[j]), 1)) # compute pitch strength S[j] = pitch_strength_one_candidate(f_erbs[int(k[j]):], nL, pc[j]) return S
def _train_pca(self, training_set): """Trains and returns a LinearMachine that is trained using PCA""" data_list = (feature for client in training_set for feature in client) data = numpy.vstack(data_list) logger.info(" -> Training Linear Machine using PCA") t = bob.learn.linear.PCATrainer() machine, eigen_values = t.train(data) if isinstance(self.pca_subspace, float): cummulated = numpy.cumsum(eigen_values) / numpy.sum(eigen_values) for index in range(len(cummulated)): if cummulated[index] > self.pca_subspace: self.pca_subspace = index break self.pca_subspace = index if self.lda_subspace is not None and self.pca_subspace <= self.lda_subspace: logger.warn(" ... Extending the PCA subspace dimension from %d to %d", self.pca_subspace, self.lda_subspace + 1) self.pca_subspace = self.lda_subspace + 1 else: logger.info(" ... Limiting PCA subspace to %d dimensions", self.pca_subspace) # limit number of pcs machine.resize(machine.shape[0], self.pca_subspace) return machine
def _interpolate_write_with_cursive(glyphs, inum, theta, noise, offset_size): stack = row_stack(glyphs) ig = _rnd_interpolate(stack, len(glyphs)*inum, ordered=True) gamma = theta + cumsum((1.0-2.0*random(len(ig)))*noise) dd = column_stack((cos(gamma), sin(gamma)))*offset_size a = ig + dd b = ig + dd[:,::-1]*array((1,-1)) return a, b
def get_frag_by_loc( cooler_file, loci, dim, zoomout_level=0, balanced=True, padding=None, percentile=100.0, ignore_diags=0, no_normalize=False ): with h5py.File(cooler_file, 'r') as f: c = get_cooler(f, zoomout_level) # Calculate the offsets once resolution = c.info['bin-size'] chromsizes = np.ceil(c.chromsizes / resolution).astype(int) offsets = np.cumsum(chromsizes) - chromsizes fragments = collect_frags( c, loci, dim, resolution, offsets, padding=padding, balanced=balanced, percentile=percentile, ignore_diags=ignore_diags, no_normalize=no_normalize ) return fragments
def abs2genomic(chromsizes, start_pos, end_pos): abs_chrom_offsets = np.r_[0, np.cumsum(chromsizes.values)] cid_lo, cid_hi = np.searchsorted(abs_chrom_offsets, [start_pos, end_pos], side='right') - 1 rel_pos_lo = start_pos - abs_chrom_offsets[cid_lo] rel_pos_hi = end_pos - abs_chrom_offsets[cid_hi] start = rel_pos_lo for cid in range(cid_lo, cid_hi): yield cid, start, chromsizes[cid] start = 0 yield cid_hi, start, rel_pos_hi
def _create_edges(self, y, order='tail'): y.sort(order=order) _docs, _counts = np.unique(y[order], return_counts=True) counts = np.zeros(self.n_docs) counts[_docs] = _counts docs = np.ascontiguousarray( np.concatenate(([0], np.cumsum(counts))), dtype=np.intc) edges = np.ascontiguousarray(y['index'].flatten(), dtype=np.intc) return docs, edges
def _compress_svd_l(self, rank, relerr, svdfunc): """Compresses the MPA in place from right to left using SVD; yields a right-canonical state See :func:`~compress` for more details and arguments. """ assert rank > 0, "Cannot compress to rank={}".format(rank) assert (relerr is None) or ((0. <= relerr) and (relerr <= 1.)), \ "relerr={} not allowed".format(relerr) for site in range(len(self) - 1, 0, -1): ltens = self._lt[site] matshape = (ltens.shape[0], -1) if relerr is None: u, sv, v = svdfunc(ltens.reshape(matshape), rank) rank_t = len(sv) else: u, sv, v = svd(ltens.reshape(matshape)) svsum = np.cumsum(sv) / np.sum(sv) rank_relerr = np.searchsorted(svsum, 1 - relerr) + 1 rank_t = min(ltens.shape[0], v.shape[0], rank, rank_relerr) yield sv, rank_t newtens = (matdot(self._lt[site - 1], u[:, :rank_t] * sv[None, :rank_t]), v[:rank_t, :].reshape((rank_t, ) + ltens.shape[1:])) self._lt.update(slice(site - 1, site + 1), newtens, canonicalization=(None, 'right')) yield np.sum(np.abs(self._lt[0])**2)
