Optimize interaction matrix scoring kernel

finches/epsilon_calculation.py

@@ -14,6 +14,38 @@
 from finches.utils import matrix_manipulation
 from finches import epsilon_stateless
 
+
+def get_attractive_repulsive_matrixes(matrix, null_interaction_baseline):
+    return epsilon_stateless.get_attractive_repulsive_matrices(matrix, null_interaction_baseline)
+
+
+def mask_matrix(matrix, column_mask):
+    return epsilon_stateless.mask_matrix(matrix, column_mask)
+
+
+def masked_matrix(matrix, column_mask):
+    return epsilon_stateless.mask_matrix(matrix, column_mask)
+
+
+def flatten_matrix_to_vector(matrix, orientation=0):
+    return np.mean(matrix, axis=orientation)
+
+
+def get_sequence_epsilon_vectors(sequence1,
+                                 sequence2,
+                                 X,
+                                 prefactor=None,
+                                 null_interaction_baseline=None,
+                                 use_charge_weighting=True,
+                                 use_aliphatic_weighting=True):
+    return epsilon_stateless.get_sequence_epsilon_vectors(sequence1,
+                                                          sequence2,
+                                                          X,
+                                                          charge_prefactor=prefactor,
+                                                          null_interaction_baseline=null_interaction_baseline,
+                                                          use_charge_weighting=use_charge_weighting,
+                                                          use_aliphatic_weighting=use_aliphatic_weighting)
+
 # -------------------------------------------------------------------------------------------------
 class InteractionMatrixConstructor:
 
@@ -691,14 +723,14 @@ def calculate_epsilon_value(self,
 
         """
 
-        # note this runs charge_prefactor and null_interaction_baseline
-        # as null which means the default values associated with this
-        # object will be used
-        return epsilon_stateless.get_sequence_epsilon_value(sequence1,
-                                          sequence2,
-                                          self,
-                                          use_charge_weighting=use_charge_weighting,
-                                          use_aliphatic_weighting=use_aliphatic_weighting)
+        w_matrix = self.calculate_weighted_pairwise_matrix(sequence1,
+                                                           sequence2,
+                                                           use_charge_weighting=use_charge_weighting,
+                                                           use_aliphatic_weighting=use_aliphatic_weighting)
+
+        transformed = w_matrix - (2.0 * self.null_interaction_baseline)
+        transformed[w_matrix == self.null_interaction_baseline] -= self.null_interaction_baseline
+        return float(np.sum(transformed) / transformed.shape[1])
 
     ## ------------------------------------------------------------------------------
     ## 
@@ -871,3 +903,5 @@ def __matrix2eps(in_matrix):
 
         return (everything, seq2_indices, seq1_indices)
 
+
+Interaction_Matrix_Constructor = InteractionMatrixConstructor

finches/forcefields/mpipi.py

@@ -600,8 +600,12 @@ def wang_frenkel(r, sigma_ij, epsilon_ij, mu_ij=2, nu_ij=1):
     main_term3 = np.power(np.power(R_ij/r, 2*mu_ij) - 1, 2*nu_ij)
 
     return main_term1*main_term2*main_term3
-
-
 
 
+def mpipi_model(version='Mpipi_GGv1', *args, **kwargs):
+    version_aliases = {
+        'mPiPi_GGv1': 'Mpipi_GGv1',
+        'mPiPi_original': 'Mpipi_original',
+    }
+    return Mpipi_model(version=version_aliases.get(version, version), *args, **kwargs)
 

finches/parsing_aminoacid_sequences.py

@@ -15,6 +15,7 @@
 """
 import numpy as np
 from finches import sequence_tools
+from finches.utils import matrix_manipulation
 
 # new characters for PIMMS aliphatic groups 
 aliphatic_group1 = {'A':'a', 'L':'l', 'M':'m', 'I':'i', 'V':'v'}
@@ -95,103 +96,8 @@ def get_charge_weighted_mask(sequence1, sequence2):
         matrix here, but this function does return
     """
 
-    #
-    # NB - this could be rewritten in Cython for improved performance...
-    #
-
-    # nb - hardcoded for now but could and probably should be altered to enable
-    # pH-dependent effects in the future
-    charges = ['R','K','E','D']
-
-    attractive_matrix = []
-    repulsive_matrix = []
-
-    n2 = len(sequence2)
-
-    # cycle through each residue
-    for i,r1 in enumerate(sequence1):
-        tmp_attractive = []
-        tmp_repulsive = []
-
-        # if r1 is charged
-        if r1 in charges:
-
-            # cycle through each residue in sequence 2
-            for j,r2 in enumerate(sequence2):
-
-                # initialize
-                w_attractive = 0
-                w_repulsive = 0
-
-                # if the second residue is charged
-                if r2 in charges: 
-
-                    # this generates a string of max 6 residues (for terminal residues 5 or 4 residues)
-                    # which is basically a concatenated fragment 
-                    l_resis = sequence_tools.get_neighbors_window_of3(i,sequence1) + sequence_tools.get_neighbors_window_of3(j,sequence2)
-
-                    # for that fragment, calculate the local fcr and ncpr
-                    [local_fcr, local_ncpr] = sequence_tools.calculate_FCR_and_NCPR(l_resis)
-
-                    # calculate the charge weight as |NCPR/FRC|. This means in one limit charge_weight goes
-                    # to 1 if the fragment is all the same type of charged residues, and goes to 0 if the
-                    # if the fragment is neutral, regardless of the fraction of charged residues.
-                    chrg_weight = np.abs(local_ncpr / local_fcr)
-
-                    w_repulsive = chrg_weight
-
-
-                    # alternative implementation - to move elsewhere at some point, but TL/DR was worse
-                    # but ALSO SLOWER! Win win!
-                    """
-
-                    frag1 = sequence_tools.get_neighbors_window_of3(i, sequence1)
-                    frag2 = sequence_tools.get_neighbors_window_of3(j, sequence2)
-
-                    [f1_fcr, f1_ncpr]  = sequence_tools.calculate_FCR_and_NCPR(frag1)
-                    [f2_fcr, f2_ncpr]  = sequence_tools.calculate_FCR_and_NCPR(frag2)
-
-                    # if both fragments contain only one type of charged residue
-                    if abs(f1_ncpr) == f1_fcr and abs(f2_ncpr) == f2_fcr:
-
-                        # calculate charge weight 
-                        q1q2 = f1_ncpr*f2_ncpr
-
-                        # opposite charge clusters
-                        if q1q2 < 0:
-
-                            # negative value (attractive) - max value = 1 (|q1| and |q2| <= 1)
-                            w_attractive = abs(q1q2)
-
-                        else:
-
-                            # positive value (repulsive) -  max value = 1 (|q1| and |q1| <= 1)
-                            w_repulsive = q1q2 
-                    """
-
-                # w_attractive and w_repulsive are 0 unless both fragements only possess the same
-                # type of charged residues
-                tmp_attractive.append(w_attractive)
-                tmp_repulsive.append(w_repulsive)
-
-
-        else:
-
-            # if r1 was not charged, create an empty vector
-            tmp_attractive = [0]*n2
-            tmp_repulsive = [0]*n2
-
-        attractive_matrix.append(tmp_attractive)
-        repulsive_matrix.append(tmp_repulsive)
-
-    # Assert matrices are the right shape
-    attractive_matrix = np.array(attractive_matrix)
-    repulsive_matrix = np.array(repulsive_matrix)
-
-    assert attractive_matrix.shape == (len(sequence1), len(sequence2))
-    assert repulsive_matrix.shape == (len(sequence1), len(sequence2))
-
-    return attractive_matrix, repulsive_matrix
+    attractive_matrix, repulsive_matrix = matrix_manipulation.charge_weighted_mask(sequence1, sequence2)
+    return np.asarray(attractive_matrix), np.asarray(repulsive_matrix)
 
 
 
@@ -281,29 +187,15 @@ def get_aliphatic_weighted_mask(sequence1, sequence2):
         returns a 2D MATRIX mask the same shape of (len(sequence1), len(sequence2))
 
     """
-    multiplier_weighting = {'1_1':1, '1_2':1, '1_3':1, '2_1':1, '3_1':1,
-                            '2_2':1.5, '2_3':1.5, '3_2':1.5, 
-                            '3_3':3}
-
-
-    ali_mask1 = get_aliphatic_groups(sequence1)
-    ali_mask2 = get_aliphatic_groups(sequence2)
-    n2 = len(sequence2)
-    matrix = []
-    for i,v1 in enumerate(ali_mask1):
-        tmp = [] 
-        if v1 > 0:
-            for j,v2 in enumerate(ali_mask2):  
-                if v2 > 0: 
-                    tmp.append(multiplier_weighting[f'{v1}_{v2}'])
-                else:
-                    tmp.append(1)
-        else:
-            tmp = [1]*n2
-
-        matrix.append(tmp)
-
-    return np.array(matrix) 
+    multiplier_weighting = np.ones((4, 4), dtype=float)
+    multiplier_weighting[2, 2] = 1.5
+    multiplier_weighting[2, 3] = 1.5
+    multiplier_weighting[3, 2] = 1.5
+    multiplier_weighting[3, 3] = 3.0
+
+    ali_mask1 = np.asarray(get_aliphatic_groups(sequence1), dtype=np.int8)
+    ali_mask2 = np.asarray(get_aliphatic_groups(sequence2), dtype=np.int8)
+    return multiplier_weighting[ali_mask1[:, None], ali_mask2[None, :]]
 
 
 ## ---------------------------------------------------------------------------

finches/tests/test_FH_diagrams.py

@@ -9,10 +9,10 @@
 from finches.forcefields.calvados import calvados_model
 from finches import epsilon_calculation
 
-from ..test_data.test_sequences import test_sequences, t0
+from .test_data.test_sequences import test_sequences, t0
 
 # test are done in the context with the mPiPi_GGv1 model
-L_model = mPiPi_model('mPiPi_GGv1')
+L_model = mpipi_model('mPiPi_GGv1')
 X_local = epsilon_calculation.Interaction_Matrix_Constructor(L_model)
 
 ############################################################################################
@@ -28,4 +28,3 @@
 #
 def build_DIELECTRIC_dependent_phase_diagrams():
     pass 
-

finches/tests/test_data/update_test_data.py

@@ -11,8 +11,9 @@
 import os
 
 import finches
+from finches import epsilon_calculation
 
-from test_sequences import test_sequences, test_condition_dict
+from .test_sequences import test_sequences, test_condition_dict, t0
 
 # ..........................................................................................
 #
@@ -78,7 +79,7 @@ def write_test_weighted_matrix(filepath, model, model_name):
         otl2.append(ot)
         print(ot.shape)
 
-    np.savez(data_file, DEFAULT=otl, NOCHARGE=otl1,  NOuse_aliphatic_weighting=otl2)
+    np.savez(data_file, DEFAULT=otl, NOCHARGE=otl1, NOALIPHATICS=otl2)
 
 #..........................................................................................
 #
@@ -101,7 +102,8 @@ def write_test_matrix_manipulation(filepath, model):
     otla, otlr = get_attractive_repulsive_matrixes(all_manipulated['test_matrix'],-0.15)
     all_manipulated['attractive_repulsive_matrixes'] = otla, otlr
 
-    mask = np.random.choice([0, 1], size=test_matrix.shape)
+    rng = np.random.default_rng(0)
+    mask = rng.choice([0, 1], size=test_matrix.shape)
     all_manipulated['bionary_mask'] = mask
 
     out_masked = epsilon_calculation.mask_matrix(test_matrix, mask)
@@ -200,4 +202,3 @@ def write_FH_out_data(filepath, model, model_name):
 
 
 
-