gistfile1.py

import sys

sys.path.append("/lhome/home/jastrzebski/mol2vec/mol2vec")
import logging
logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.DEBUG)
from sklearn.cross_validation import StratifiedKFold
from experiments.utils import wac_score
from sklearn.svm import SVC
from sklearn.grid_search import GridSearchCV
from training_data.datasets import *
from sklearn.preprocessing import MaxAbsScaler
from collections import defaultdict

def _calculate_jaccard_kernel(X1T, X2T):
    X1T_sums = np.array(X1T.sum(axis=1))
    X2T_sums = np.array(X2T.sum(axis=1))
    K = X1T.dot(X2T.T)

    if hasattr(K, "toarray"):
        K = K.toarray()

    K2 = -(K.copy())
    K2 += (X1T_sums.reshape(-1, 1))
    K2 += (X2T_sums.reshape(1, -1))
    K = K / K2

    return K


# Construct triples
active_file = "klek/5-HT1a_actives_0_KRFP_coord.dat"
inactive_file = "klek/5-HT1a_inactives_0_KRFP_coord.dat"
data = []
indptr = [0]
col_ptr = []
id = 0
y = []
KLEK_SIZE = 4860  # Krzysiek uses 4860 bits of KR

for f in [active_file, inactive_file]:
    label = int("_actives" in f)
    examples = 0
    for line in open(f, "r"):
        examples += 1
        tokens = line.split()
        for token in tokens[1:]:  # Skip CHEMBL id
            feature_name, feature_val = token.split(":")
            v1, v2 = feature_name.split(",")
            v1, v2 = int(v1), int(v2)
            assert v1 < KLEK_SIZE and v2 < KLEK_SIZE, "Indexing must be unique"
            id += 1
            data.append(int(feature_val))
            col_ptr.append(v1 * KLEK_SIZE + v2)

        indptr.append(id)  # indptr[i+1] - index of last elemnt in ith column
    y += [label] * examples

data = np.array(data).astype("int")
X = scipy.sparse.csr_matrix((data, col_ptr, indptr))

C_min = -5
C_max = 9
rng = 777
max_iter = 6 * 10 ** 5
n_folds = 5
X_abs = MaxAbsScaler().fit_transform(X)
K = _calculate_jaccard_kernel(X_abs, X_abs)
K[np.isnan(K)] = 0.

param_grid = {"C": [10 ** i for i in range(C_min, C_max)]}
estimator = SVC(kernel="precomputed", class_weight="balanced", random_state=rng, max_iter=max_iter)
cv = GridSearchCV(estimator=estimator, n_jobs=1,
                  cv=StratifiedKFold(y, random_state=rng, shuffle=True, n_folds=n_folds),
                  param_grid=param_grid, scoring=wac_scoring)
cv.fit(K, y)

logger.info(cv.best_score_)
logger.info(cv.grid_scores_)