From 9c0c8a9d910913b3952353806eb82b367fe0a48a Mon Sep 17 00:00:00 2001 From: xwanaf <66076868+xwanaf@users.noreply.github.com> Date: Wed, 12 Jul 2023 16:31:02 +0800 Subject: [PATCH] Update utils_pyRCTD.py --- src/utils_pyRCTD.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/utils_pyRCTD.py b/src/utils_pyRCTD.py index 276f13d..01df419 100644 --- a/src/utils_pyRCTD.py +++ b/src/utils_pyRCTD.py @@ -157,7 +157,7 @@ def psd(H): # P = eig[1] @ np.diag(np.clip(eig[0], a_min = epsilon, a_max = eig[0].max() + 10)) @ eig[1].T return P -def solveWLS(S,B,S_mat,initialSol, nUMI, bulk_mode = False, constrain = False, likelihood_vars = None): +def solveWLS(S,B,S_mat,initialSol, nUMI, bulk_mode = False, constrain = False, likelihood_vars = None, solver = 'osqp'): solution = initialSol.copy() solution[solution < 0] = 0 prediction = np.absolute(S @ solution) @@ -176,13 +176,13 @@ def solveWLS(S,B,S_mat,initialSol, nUMI, bulk_mode = False, constrain = False, l bzero = -solution alpha = 0.3 if constrain: - solution = solution + alpha*solve_qp(np.array(D_mat),-np.array(d_vec),-np.array(A),-np.array(bzero), np.ones(solution.shape[0]), 1 - solution.sum(), solver="osqp") + solution = solution + alpha*solve_qp(np.array(D_mat),-np.array(d_vec),-np.array(A),-np.array(bzero), np.ones(solution.shape[0]), 1 - solution.sum(), solver=solver) else: - solution = solution + alpha*solve_qp(np.array(D_mat),-np.array(d_vec),-np.array(A),-np.array(bzero), solver="osqp") + solution = solution + alpha*solve_qp(np.array(D_mat),-np.array(d_vec),-np.array(A),-np.array(bzero), solver=solver) return solution def solveIRWLS_weights(S,B,nUMI, OLS=False, constrain = True, verbose = False, - n_iter = 50, MIN_CHANGE = .001, bulk_mode = False, solution = None, loggings = None, likelihood_vars = None): + n_iter = 50, MIN_CHANGE = .001, bulk_mode = False, solution = None, loggings = None, likelihood_vars = None, solver = 'osqp'): if not bulk_mode: K_val = likelihood_vars['K_val'] B = np.copy(B) @@ -193,7 +193,7 @@ def solveIRWLS_weights(S,B,nUMI, OLS=False, constrain = True, verbose = False, change = 1 changes = [] while (change > MIN_CHANGE) & (iterations < n_iter): - new_solution = solveWLS(S,B,S_mat,solution, nUMI,constrain=constrain, bulk_mode = bulk_mode, likelihood_vars = likelihood_vars) + new_solution = solveWLS(S,B,S_mat,solution, nUMI,constrain=constrain, bulk_mode = bulk_mode, likelihood_vars = likelihood_vars, solver = solver) change = np.linalg.norm(new_solution-solution, 1) if verbose: loggings.info('Change: {}'.format(change)) @@ -208,8 +208,12 @@ def decompose_full_ray(args): bulk_mode = False verbose = False n_iter = 50 - results = solveIRWLS_weights(cell_type_profiles,bead,nUMI,OLS = OLS, constrain = constrain, - verbose = verbose, n_iter = n_iter, MIN_CHANGE = MIN_CHANGE, bulk_mode = bulk_mode, loggings = loggings, likelihood_vars = likelihood_vars) + try: + results = solveIRWLS_weights(cell_type_profiles,bead,nUMI,OLS = OLS, constrain = constrain, + verbose = verbose, n_iter = n_iter, MIN_CHANGE = MIN_CHANGE, bulk_mode = bulk_mode, loggings = loggings, likelihood_vars = likelihood_vars, solver = 'osqp') + except: + results = solveIRWLS_weights(cell_type_profiles,bead,nUMI,OLS = OLS, constrain = constrain, + verbose = verbose, n_iter = n_iter, MIN_CHANGE = MIN_CHANGE, bulk_mode = bulk_mode, loggings = loggings, likelihood_vars = likelihood_vars, solver = 'cvxopt') return results @@ -690,4 +694,4 @@ def check_pairs_type(cell_type_profiles, bead, UMI_tot, score_mat, min_score, my if all_pairs_class and not all_pairs and (len(other_class) > 1): for ty in other_class[2:len(other_class)]: singlet_score = min(singlet_score, get_singlet_score(cell_type_profiles, bead, UMI_tot, ty, constrain, MIN_CHANGE = MIN_CHANGE, loggings = loggings, likelihood_vars = likelihood_vars)) - return {'all_pairs': all_pairs, 'all_pairs_class': all_pairs_class, 'singlet_score': singlet_score} \ No newline at end of file + return {'all_pairs': all_pairs, 'all_pairs_class': all_pairs_class, 'singlet_score': singlet_score}