The Maximisation step of the CEM algorithm
435452
def maximise(self, **kwargs):
        
        self.scorer.write_partition(self.partition)
        self.scorer.analyse_cache_dir(**kwargs)
        self.likelihood = self.scorer.get_partition_score(self.partition)
        self.scorer.clean_cache()
        changed = self.get_changed(self.partition, self.prev_partition)
        self.update_perlocus_likelihood_objects(self.partition, changed)
        return self.partition, self.likelihood, sum(inst.get_likelihood() for inst in self.insts)
585804
private void mmuli(int n, double[] A, double[] B, double[] C) {
        for (int i = 0; i < n * n; i++) {
            C[i] = 0;
        }

        for (int j = 0; j < n; j++) {
            int jn = j * n;
            for (int k = 0; k < n; k++) {
                int kn = k * n;
                double bkjn = B[k + jn];
                for (int i = 0; i < n; i++) {
                    C[i + jn] += A[i + kn] * bkjn;
                }
            }
        }
    }
508264
def cv_precompute(self, mask, b):
        '''
        

        '''

        # Get current chunk and mask outliers
        m1 = self.get_masked_chunk(b)
        flux = self.fraw[m1]
        K = GetCovariance(self.kernel, self.kernel_params,
                          self.time[m1], self.fraw_err[m1])
        med = np.nanmedian(flux)

        # Now mask the validation set
        M = lambda x, axis = 0: np.delete(x, mask, axis=axis)
        m2 = M(m1)
        mK = M(M(K, axis=0), axis=1)
        f = M(flux) - med

        # Pre-compute the matrices
        A = [None for i in range(self.pld_order)]
        B = [None for i in range(self.pld_order)]
        for n in range(self.pld_order):
            # Only compute up to the current PLD order
            if self.lam_idx >= n:
                X2 = self.X(n, m2)
                X1 = self.X(n, m1)
                A[n] = np.dot(X2, X2.T)
                B[n] = np.dot(X1, X2.T)
                del X1, X2

        if self.transit_model is None:
            C = 0
        else:
            C = np.zeros((len(m2), len(m2)))
            mean_transit_model = med * \
                np.sum([tm.depth * tm(self.time[m2])
                        for tm in self.transit_model], axis=0)
            f -= mean_transit_model
            for tm in self.transit_model:
                X2 = tm(self.time[m2]).reshape(-1, 1)
                C += tm.var_depth * np.dot(X2, X2.T)
                del X2

        return A, B, C, mK, f, m1, m2
440714
function () {
			var k, ent,
			hash = {
				max: { x: -Infinity, y: -Infinity },
				min: { x: Infinity, y: Infinity }
			},
			coords = {
				max: { x: -Infinity, y: -Infinity },
				min: { x: Infinity, y: Infinity }
			};

      //Using broad phase hash to speed up the computation of boundaries.
			for (var h in this.map) {
				if (!this.map[h].length) continue;

        //broad phase coordinate
				var i= h>>16,
					j=(h<<16)>>16;
				if (j<0) { i = i^-1 }
				if (i >= hash.max.x) {
					hash.max.x = i;
					for (k in this.map[h]) {
						ent = this.map[h][k];
						//make sure that this is a Crafty entity
						if (typeof ent == 'object' && 'requires' in ent) {
							coords.max.x = Math.max(coords.max.x, ent.x + ent.w);
						}
					}
				}
				if (i <= hash.min.x) {
					hash.min.x = i;
					for (k in this.map[h]) {
						ent = this.map[h][k];
						if (typeof ent == 'object' && 'requires' in ent) {
							coords.min.x = Math.min(coords.min.x, ent.x);
						}
					}
				}
				if (j >= hash.max.y) {
					hash.max.y = j;
					for (k in this.map[h]) {
						ent = this.map[h][k];
						if (typeof ent == 'object' && 'requires' in ent) {
							coords.max.y = Math.max(coords.max.y, ent.y + ent.h);
						}
					}
				}
				if (j <= hash.min.y) {
					hash.min.y = j;
					for (k in this.map[h]) {
						ent = this.map[h][k];
						if (typeof ent == 'object' && 'requires' in ent) {
							coords.min.y = Math.min(coords.min.y, ent.y);
						}
					}
				}
			}

			return coords;
		}
457325
def ccmod_xstep(k):
    

    YU = mp_D_Y - mp_D_U[k]
    b = mp_ZSf[k] + mp_drho * sl.rfftn(YU, None, mp_cri.axisN)
    Xf = sl.solvedbi_sm(mp_Zf[k], mp_drho, b, axis=mp_cri.axisM)
    mp_D_X[k] = sl.irfftn(Xf, mp_cri.Nv, mp_cri.axisN)
1153674
def _maximization_step(X, posteriors):
    

    ### Prior probabilities or class weights
    sum_post_proba = np.sum(posteriors, axis=0)
    prior_proba = sum_post_proba / (sum_post_proba.sum() + Epsilon)
    
    ### means
    means = np.dot(posteriors.T, X) / (sum_post_proba[:, np.newaxis] + Epsilon)

    ### covariance matrices
    n_components = posteriors.shape[1]
    n_features = X.shape[1]
    covars = np.empty(shape=(n_components, n_features, n_features), dtype=float)
  
    for i in range(n_components):
        post_i = posteriors[:, i]
        mean_i = means[i]
        diff_i = X - mean_i

        with np.errstate(under='ignore'):
            covar_i = np.dot(post_i * diff_i.T, diff_i) / (post_i.sum() + Epsilon)
        covars[i] = covar_i + Lambda * np.eye(n_features)


    _validate_params(prior_proba, means, covars)
    return(prior_proba, means, covars)
1165732
def _preoptimize_model(self):
        
        gaussian_model = LLT(self.data, integ=self.integ, target=self.target)
        gaussian_model.fit()
        self.latent_variables.z_list[0].start = gaussian_model.latent_variables.get_z_values()[1]
        self.latent_variables.z_list[1].start = gaussian_model.latent_variables.get_z_values()[2]

        if self.model_name2 == 't':

            def temp_function(params):
                return -np.sum(ss.t.logpdf(x=self.data, df=np.exp(params[0]), 
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2])))

            p = optimize.minimize(temp_function,np.array([2.0, 0.0, -1.0]),method='L-BFGS-B')
            self.latent_variables.z_list[2].start = p.x[2]
            self.latent_variables.z_list[3].start = p.x[0]

        elif self.model_name2 == 'Skewt':

            def temp_function(params):
                return -np.sum(fam.Skewt.logpdf_internal(x=self.data,df=np.exp(params[0]),
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2]),gamma=np.exp(params[3])))

            p = optimize.minimize(temp_function,np.array([2.0, 0.0, -1.0, 0.0]),method='L-BFGS-B')
            self.latent_variables.z_list[2].start = p.x[3]
            self.latent_variables.z_list[3].start = p.x[2]
            self.latent_variables.z_list[4].start = p.x[0]

        return gaussian_model.latent_variables
384844
def fit(self, X):
        
        self.centers_, self.labels_, self.sse_arr_, self.n_iter_ = \
              _kmeans(X, self.n_clusters, self.max_iter, self.n_trials, self.tol)
304400
def _pso(self, n_particles, n_iterations, optimizer):

        

        pso = ParticleSwarmOptimizer(optimizer, low=self._lower_limit, high=self._upper_limit, particleCount=n_particles)

        gBests = pso._optimize(maxIter=n_iterations,standard_dev=self._pso_convergence_standardDEV)

        likelihoods = [particle.fitness for particle in gBests]
        ind = np.argmax(likelihoods)

        return gBests[ind].position
461514
def partition_inference(tctc_mat, comp, tau, sigma, kappa):
     
    """
    communityinfo = {}
    communityinfo['community'] = []
    communityinfo['start'] = np.empty(0)
    communityinfo['end'] = np.empty(0)
    communityinfo['size'] = np.empty(0)
    for i, tcomp in enumerate(comp):
        # This can go in parallel loop
        if len(tcomp) > 0:
            for traj in tcomp:
                # Check it does not already exist.
                ignore = 0
                preexisting = 0
                if i != 0:
                    cutoff = i-1-kappa
                    if cutoff < 0:
                        cutoff = 0
                    
                    if np.any(np.sum(np.sum(tctc_mat[traj, :, cutoff:i][:, traj], axis=0), axis=0) == np.power(len(traj), 2)):
                        # Make sure that a small trajectory could exist 
                        for checknode in np.where(communityinfo['end']>=cutoff)[0]:
                            if traj == communityinfo['community'][checknode]:
                                ignore = 1
                        if ignore == 0: 
                            for checknode in np.where(communityinfo['end']>=cutoff)[0]:
                                if set(communityinfo['community'][checknode]).issuperset(traj):
                                    preexisting = 1                              
                if ignore == 0:
                    # Check how long it continues
                    # For efficiency, increase in blocks
                    approxmaxlength = tau*2
                    a = np.sum(
                        np.sum(tctc_mat[traj, :, i:i+approxmaxlength][:, traj], axis=0), axis=0)
                    if len(traj)*len(traj)*approxmaxlength == a.sum():
                        ok = 0
                        ii = 1
                        while ok == 0:
                            b = np.sum(np.sum(
                                tctc_mat[traj, :, i+(approxmaxlength*ii):i+(approxmaxlength*(ii+1))][:, traj], axis=0), axis=0)
                            a = np.append(a, b)
                            if len(traj)*len(traj)*approxmaxlength != b.sum():
                                ok = 1
                            if i+(approxmaxlength*(ii+1)) > tctc_mat.shape[-1]:
                                ok = 1
                            ii += 1
                    a = np.where(a == np.power(len(traj), 2))[0]
                    # Add an additional value that is false in case end of time series
                    if len(a) == 1:
                        stopind = i + 1
                    else:
                        a = np.append(a, a.max()+kappa+2)
                        # Find the stop index (if stopind = 4 and start = 0, then tctc_mat[:,:,start:stopind]==1)
                        stopind = i + np.split(a, np.where(
                                np.diff(a) > kappa+1)[0]+1)[0][-1] + 1
                    # Add trajectory to dictionary
                    if ((stopind - i) >= tau or preexisting == 1) and len(traj) >= sigma:
                        communityinfo['community'].append(sorted(traj))
                        communityinfo['start'] = np.append(communityinfo['start'], int(i))
                        communityinfo['end'] = np.append(
                            communityinfo['end'], int(stopind))
                        communityinfo['size'] = np.append(communityinfo['size'], len(traj))

    communityinfo = pd.DataFrame(communityinfo)

    communityinfo['start'] = communityinfo['start'].astype(int)
    communityinfo['end'] = communityinfo['end'].astype(int)
    # First check that there is not already a trajectory that is ongoing
    badrows = []
    for v in communityinfo.iterrows():
        skipselrule = (communityinfo['end'] == v[1]['end'])
        for u in communityinfo[skipselrule].iterrows():
            a = 1
            if u[1]['start'] > v[1]['start'] and sorted(u[1]['community']) == sorted(v[1]['community']):
                badrows.append(u[0])
    communityinfo = communityinfo.drop(badrows)

    # Then see if any subset trajectory can be placed earlier in time.
    for v in communityinfo.iterrows():
        skipselrule = (communityinfo['end'] <= v[1]['start']) & (
            communityinfo['end']+kappa >= v[1]['start'])
        for u in communityinfo[skipselrule].iterrows():
            a = 1
            if set(u[1]['community']).issuperset(v[1]['community']):
                communityinfo.loc[v[0], 'start'] = u[1]['start']

    # It is possible to make the condition below effective_length
    communityinfo['length'] = np.array(communityinfo['end']) - np.array(communityinfo['start'])
    communityinfo = communityinfo[communityinfo['length'] >= tau]
    communityinfo = communityinfo[communityinfo['size'] >= sigma]

    # Make sure that the traj is not completely enguled by another
    badrows = []
    if kappa > 0:
        for v in communityinfo.iterrows():
            skipselrule = (communityinfo['end'] == v[1]['end']) & (
                communityinfo['start'] < v[1]['start'])
            for u in communityinfo[skipselrule].iterrows():
                if set(v[1]['community']).issubset(u[1]['community']):
                    badrows.append(v[0])
        communityinfo = communityinfo.drop(badrows)

    return communityinfo
1162081
def create(max_kl, cg_iters, line_search_iters, cg_damping, entropy_coef, vf_iters, discount_factor,
           gae_lambda=1.0, improvement_acceptance_ratio=0.1, max_grad_norm=0.5):
    
    return TrpoPolicyGradient(
        max_kl, int(cg_iters), int(line_search_iters), cg_damping, entropy_coef, vf_iters,
        discount_factor=discount_factor,
        gae_lambda=gae_lambda,
        improvement_acceptance_ratio=improvement_acceptance_ratio,
        max_grad_norm=max_grad_norm
    )
205233
def solve_mle_rev(C, tol=1e-10, maxiter=100, show_progress=False, full_output=False,
                  return_statdist=True, **kwargs):
    """"""
    M = C.shape[0]

    """Initial guess for primal-point"""
    z0 = np.zeros(2*M)
    z0[0:M] = 1.0

    """Inequality constraints"""
    # G = np.zeros((M, 2*M))
    # G[np.arange(M), np.arange(M)] = -1.0
    G = -1.0*scipy.sparse.eye(M, n=2*M, k=0)
    h = np.zeros(M)

    """Equality constraints"""
    A = np.zeros((1, 2*M))
    A[0, M] = 1.0
    b = np.array([0.0])

    """Scaling"""
    c0 = C.max()
    C = C/c0

    """Symmetric part"""
    Cs = C + C.T

    """Column sum"""
    c = C.sum(axis=0)

    if scipy.sparse.issparse(C):
        Cs = Cs.tocsr()
        c = c.A1
        A = scipy.sparse.csr_matrix(A)
        F = objective_sparse.F
        DF = objective_sparse.DFsym
        convert_solution = objective_sparse.convert_solution
    else:
        F = objective_dense.F
        DF = objective_dense.DF
        convert_solution = objective_dense.convert_solution

    """PDIP iteration"""
    res = primal_dual_solve(F, z0, DF, A, b, G, h,
                            args=(Cs, c),
                            maxiter=maxiter, tol=tol,
                            show_progress=show_progress,
                            full_output=full_output)
    if full_output:
        z, info = res
    else:
        z = res
    pi, P = convert_solution(z, Cs)

    result = [P]
    if return_statdist:
        result.append(pi)
    if full_output:
        result.append(info)

    return tuple(result) if len(result) > 1 else result[0]
282040
def _ces_distance_emd(unique_C1, unique_C2):
    
    # Get the pairwise distances between the concepts in the unpartitioned and
    # partitioned CESs.
    distances = np.array([
        [concept_distance(i, j) for j in unique_C2] for i in unique_C1
    ])
    # We need distances from all concepts---in both the unpartitioned and
    # partitioned CESs---to the null concept, because:
    # - often a concept in the unpartitioned CES is destroyed by a
    #   cut (and needs to be moved to the null concept); and
    # - in certain cases, the partitioned system will have *greater* sum of
    #   small-phi, even though it has less big-phi, which means that some
    #   partitioned-CES concepts will be moved to the null concept.
    distances_to_null = np.array([
        concept_distance(c, c.subsystem.null_concept)
        for ces in (unique_C1, unique_C2) for c in ces
    ])
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Now we make the distance matrix, which will look like this:
    #
    #        C1       C2     0
    #    +~~~~~~~~+~~~~~~~~+~~~+
    #    |        |        |   |
    # C1 |   X    |    D   |   |
    #    |        |        |   |
    #    +~~~~~~~~+~~~~~~~~+ D |
    #    |        |        | n |
    # C2 |   D'   |    X   |   |
    #    |        |        |   |
    #    +~~~~~~~~+~~~~~~~~+~~~|
    #  0 |        Dn'      | X |
    #    +~~~~~~~~~~~~~~~~~~~~~+
    #
    # The diagonal blocks marked with an X are set to a value larger than any
    # pairwise distance between concepts. This ensures that concepts are never
    # moved to another concept within their own CES; they must always go either
    # from one CES to another, or to the null concept N. The D block is filled
    # with the pairwise distances between the two CESs, and Dn is filled with
    # the distances from each concept to the null concept.
    N, M = len(unique_C1), len(unique_C2)
    # Add one to the side length for the null concept distances.
    distance_matrix = np.empty([N + M + 1] * 2)
    # Ensure that concepts are never moved within their own CES.
    distance_matrix[:] = np.max(distances) + 1
    # Set the top-right block to the pairwise CES distances.
    distance_matrix[:N, N:-1] = distances
    # Set the bottom-left block to the same, but transposed.
    distance_matrix[N:-1, :N] = distances.T
    # Do the same for the distances to the null concept.
    distance_matrix[-1, :-1] = distances_to_null
    distance_matrix[:-1, -1] = distances_to_null.T
    distance_matrix[-1, -1] = 0
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Construct the two phi distributions, with an entry at the end for the
    # null concept.
    d1 = [c.phi for c in unique_C1] + [0] * M + [0]
    d2 = [0] * N + [c.phi for c in unique_C2] + [0]
    # Calculate how much phi disappeared and assign it to the null concept.
    d2[-1] = sum(d1) - sum(d2)
    # The sum of the two signatures should be the same.
    assert utils.eq(sum(d1), sum(d2))
    # Calculate!
    return emd(np.array(d1), np.array(d2), distance_matrix)
681440
def precompute_optimzation_Y(laplacian_matrix, n_samples, relaxation_kwds):
    """"""
    relaxation_kwds.setdefault('presave',False)
    relaxation_kwds.setdefault('presave_name','pre_comp_current.npy')
    relaxation_kwds.setdefault('verbose',False)
    if relaxation_kwds['verbose']:
        print ('Making Lk and nbhds')
    Lk_tensor, nbk, si_map = \
        compute_Lk(laplacian_matrix, n_samples, relaxation_kwds['subset'])
    if relaxation_kwds['presave']:
        raise NotImplementedError('Not yet implemented presave')
    return { 'Lk': Lk_tensor, 'nbk': nbk, 'si_map': si_map }
1084635
def adapt(self, d, x):
        
        y = np.dot(self.w, x)
        e = d - y
        R1 = np.dot(np.dot(np.dot(self.R,x),x.T),self.R)
        R2 = self.mu + np.dot(np.dot(x,self.R),x.T)
        self.R = 1/self.mu * (self.R - R1/R2)
        dw = np.dot(self.R, x.T) * e
        self.w += dw
1150485
def adapt(self, d, x):
        
        # create input matrix and target vector
        self.x_mem[:,1:] = self.x_mem[:,:-1]
        self.x_mem[:,0] = x
        self.d_mem[1:] = self.d_mem[:-1]
        self.d_mem[0] = d
        # estimate output and error
        self.y_mem = np.dot(self.x_mem.T, self.w)
        self.e_mem = self.d_mem - self.y_mem
        # update
        dw_part1 = np.dot(self.x_mem.T, self.x_mem) + self.ide_eps
        dw_part2 = np.linalg.solve(dw_part1, self.ide)
        dw = np.dot(self.x_mem, np.dot(dw_part2, self.e_mem))
        self.w += self.mu * dw
697161
private int getMaxOccurence(double massTo, int element_pos, int[] matrix, List<IIsotope> isoToCond_new) {
        double massIn = isoToCond_new.get(element_pos).getExactMass();
        double massToM = massTo;
        for (int i = 0; i < matrix.length; i++)
            if (i != element_pos) if (matrix[i] != 0) massToM -= isoToCond_new.get(i).getExactMass() * matrix[i];

        int value = (int) ((massToM + 1) / massIn);
        return value;
    }