<jupyter_start><jupyter_text>Get validation data<jupyter_code>data_manager = data_loader.AMPDataManager(
    '../data/unlabelled_positive.csv',
    '../data/unlabelled_negative.csv',
    min_len=MIN_LENGTH,
    max_len=MAX_LENGTH)

amp_x, amp_y = data_manager.get_merged_data()
amp_x_train, amp_x_test, amp_y_train, amp_y_test = train_test_split(amp_x, amp_y, test_size=0.1, random_state=36)
amp_x_train, amp_x_val, amp_y_train, amp_y_val = train_test_split(amp_x_train, amp_y_train, test_size=0.2, random_state=36)
# Restrict the length
ecoli_df = pd.read_csv('../data/mic_data.csv')
mask = (ecoli_df['sequence'].str.len() <= MAX_LENGTH) & (ecoli_df['sequence'].str.len() >= MIN_LENGTH)
ecoli_df = ecoli_df.loc[mask]
mic_x = pad(to_one_hot(ecoli_df['sequence']))
mic_y = ecoli_df.value

mic_x_train, mic_x_test, mic_y_train, mic_y_test = train_test_split(mic_x, mic_y, test_size=0.1, random_state=36)
mic_x_train, mic_x_val, mic_y_train, mic_y_val = train_test_split(mic_x_train, mic_y_train, test_size=0.2, random_state=36)
pos = np.vstack([amp_x_val[amp_y_val == 1], mic_x_val[mic_y_val < 1.5]])
neg = np.vstack([amp_x_val[amp_y_val == 0], mic_x_val[mic_y_val > 1.5]])
neg.shape, pos.shape
pos_amp = amp_classifier_model.predict(pos, verbose=1).reshape(len(pos))
neg_mic = mic_classifier_model.predict(neg, verbose=1).reshape(len(neg))
neg_amp = amp_classifier_model.predict(neg, verbose=1).reshape(len(neg))
pos_mic = mic_classifier_model.predict(pos, verbose=1).reshape(len(pos))
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(12, 4), sharex=True, sharey=True) 
ax1.hist(pos_amp)
ax1.set_ylabel('AMP')
ax1.set_title('Positives')
ax2.hist(neg_amp)
ax2.set_title('Negatives')

ax3.hist(pos_mic)
ax3.set_ylabel('MIC')
ax4.hist(neg_mic)
plt.show()
pos = np.vstack([pos] * 64).reshape(-1, 25)
pos_amp = np.vstack([pos_amp] * 64).reshape(-1, 1)
pos_mic = np.vstack([pos_mic] * 64).reshape(-1, 1)
neg = np.vstack([neg] * 64).reshape(-1, 25)
neg_amp = np.vstack([neg_amp] * 64).reshape(-1, 1)
neg_mic = np.vstack([neg_mic] * 64).reshape(-1, 1)
def improve(x, model, epoch, mode):
    if mode == 'pos':
        amp = pos_amp
        mic = pos_mic
    else:
        amp = neg_mic
        mic = neg_mic
        
    encoded = encoder_model.predict(x, batch_size=5000)
    conditioned = np.hstack([
        encoded,
        np.ones((len(x), 1)),
        np.ones((len(x), 1)),
    ])
    decoded = decoder_model.predict(conditioned, batch_size=5000)
    new_peptides = np.argmax(decoded, axis=2)
    new_amp = amp_classifier_model.predict(new_peptides, batch_size=5000)
    new_mic = mic_classifier_model.predict(new_peptides, batch_size=5000)                                                                             
    
    # RELATIVE
    rel_better = new_amp > amp.reshape(-1, 1)
    rel_better = rel_better & (new_mic > mic.reshape(-1, 1))
    rel_better = np.logical_or.reduce(rel_better, axis=1)
    rel_improved = new_peptides[np.where(rel_better), :].reshape(-1, 25)
    before_rel_improve = x[np.where(rel_better), :].reshape(-1, 25)
    
    # ABSOLUTE
    abs_better = new_amp >= 0.8
    abs_better = abs_better & (new_mic > 0.5)
    abs_better = np.logical_or.reduce(abs_better, axis=1)
    abs_improved = new_peptides[np.where(abs_better), :].reshape(-1, 25)
    before_abs_improve = x[np.where(abs_better), :].reshape(-1, 25)
    
    return {
        'new_peptides': new_peptides,
        'rel_improved': rel_improved,
        'abs_improved': abs_improved,
        'before_rel_improve': before_rel_improve,
        'before_abs_improve': before_abs_improve, 
        'new_amp': new_amp,
        'new_mic': new_mic,
        }<jupyter_output><empty_output><jupyter_text>HydrAMP improve<jupyter_code>from keras.models import Model
model = models[0]
current_model_pos = {epoch: [] for epoch in range(40)} 
current_model_neg = {epoch: [] for epoch in range(40)}    
for epoch in tqdm(range(40)):
        AMPMaster = bms.load_model(f'../models/{model}/{epoch}')
        encoder_model =  AMPMaster.encoder(input_to_encoder)
        decoder_model = AMPMaster.decoder(input_to_decoder)         
        current_model_pos[epoch] = improve(pos, model, epoch, 'pos')
        current_model_neg[epoch] = improve(neg, model, epoch, 'neg')
dump(current_model_pos, f'../results/improvement_PosVal_{model}.joblib') 
dump(current_model_neg, f'../results/improvement_NegVal_{model}.joblib')<jupyter_output>0%|          | 0/40 [00:00<?, ?it/s]<jupyter_text>PepCVAE improve<jupyter_code>from keras.models import Model
model = models[1]
current_model_pos = {epoch: [] for epoch in range(40)} 
current_model_neg = {epoch: [] for epoch in range(40)}    
for epoch in tqdm(range(40)):
        AMPMaster = bms.load_model(f'../models/{model}/{epoch}')
        encoder_model =  AMPMaster.encoder(input_to_encoder)
        decoder_model = AMPMaster.decoder(input_to_decoder)
        new_act = layers.TimeDistributed(
            layers.Activation(activations.softmax),
            name='decoder_time_distribute_activation')
        decoder_model.layers.pop()
        x = new_act(decoder_model.layers[-1].output)
        decoder_model = Model(input=decoder_model.input, output=[x])            
        current_model_pos[epoch] = improve(pos, model, epoch, 'pos')
        current_model_neg[epoch] = improve(neg, model, epoch, 'neg')
dump(current_model_pos, f'../results/improvement_PosVal_{model}.joblib') 
dump(current_model_neg, f'../results/improvement_NegVal_{model}.joblib')<jupyter_output>0%|          | 0/40 [00:00<?, ?it/s]<jupyter_text>Basic improvement<jupyter_code>from keras.models import Model
model = models[2]
current_model_pos = {epoch: [] for epoch in range(40)} 
current_model_neg = {epoch: [] for epoch in range(40)}    
for epoch in tqdm(range(40)):
        AMPMaster = bms.load_model(f'../models/{model}/{epoch}')
        encoder_model =  AMPMaster.encoder(input_to_encoder)
        decoder_model = AMPMaster.decoder(input_to_decoder)
        new_act = layers.TimeDistributed(
            layers.Activation(activations.softmax),
            name='decoder_time_distribute_activation')
        decoder_model.layers.pop()
        x = new_act(decoder_model.layers[-1].output)
        decoder_model = Model(input=decoder_model.input, output=[x])            
        current_model_pos[epoch] = improve(pos, model, epoch, 'pos')
        current_model_neg[epoch] = improve(neg, model, epoch, 'neg')
dump(current_model_pos, f'../results/improvement_PosVal_{model}.joblib') 
dump(current_model_neg, f'../results/improvement_NegVal_{model}.joblib')<jupyter_output>0%|          | 0/40 [00:00<?, ?it/s]<jupyter_text>Collect results<jupyter_code>pos_final_results =  {model: {epoch: 
                          {'absolute improvement':0,
                           'relative improvement':0,
                          } for epoch in range(40)} for model in models}
neg_final_results = {model: {epoch: 
                          {'absolute improvement':0,
                           'relative improvement':0,
                          } for epoch in range(40)} for model in models}
for model in models:
    if model in ['PepCVAE', 'Basic']:
        model_results = load(f'../results/improvement_PosVal_{model}.joblib')
    else:
        model_results = load(f'../results/improvement_PosVal_{model}.joblib')
    for epoch in range(40):
        pos_final_results[model][epoch]['relative improvement'] = np.unique(
            model_results[epoch]['rel_improved'], axis=0).shape[0]
        pos_final_results[model][epoch]['absolute improvement'] = np.unique(
             model_results[epoch]['abs_improved'], axis=0).shape[0]
        pos_final_results[model][epoch]['before relative improvement'] = np.unique(
            model_results[epoch]['before_rel_improve'], axis=0).shape[0]
        pos_final_results[model][epoch]['before absolute improvement'] = np.unique(
            model_results[epoch]['before_abs_improve'], axis=0).shape[0]
        
        
for model in models:
    if model in ['PepCVAE', 'Basic']:
        model_results = load(f'../results/improvement_NegVal_{model}.joblib')
    else:
        model_results = load(f'../results/improvement_NegVal_{model}.joblib')
    for epoch in range(40):
            neg_final_results[model][epoch]['relative improvement'] = np.unique(
                model_results[epoch]['rel_improved'], axis=0).shape[0]
            neg_final_results[model][epoch]['absolute improvement'] = np.unique(
                model_results[epoch]['abs_improved'], axis=0).shape[0]
            neg_final_results[model][epoch]['before relative improvement'] = np.unique(
                model_results[epoch]['before_rel_improve'], axis=0).shape[0]
            neg_final_results[model][epoch]['before absolute improvement'] = np.unique(
                model_results[epoch]['before_abs_improve'], axis=0).shape[0]
hydra_metrics = pd.read_csv('../models/HydrAMP/metrics.csv')
pepcvae_metrics = pd.read_csv('../models/PepCVAE/metrics.csv')
basic_metrics = pd.read_csv('../models/Basic/metrics.csv')

plt.title('Relative improved')
plt.plot([pos_final_results[models[0]][epoch]['relative improvement']  for epoch in range(10, 40)], c='red', label='HydrAMP')
plt.plot([pos_final_results[models[1]][epoch]['relative improvement'] for epoch in range(10, 40)], c='orange', label='PepCVAE')
plt.plot([pos_final_results[models[2]][epoch]['relative improvement'] for epoch in range(10, 40)], c='blue', label='Basic')

plt.legend(bbox_to_anchor=(1.1, 0.5))
plt.show()
plt.title('Absolute improved')
plt.plot([pos_final_results[models[0]][epoch]['absolute improvement']  for epoch in range(10, 40)], c='red', label='HydrAMP')
plt.plot([pos_final_results[models[1]][epoch]['absolute improvement'] for epoch in range(10, 40)], c='orange', label='PepCVAE')
plt.plot([pos_final_results[models[2]][epoch]['absolute improvement'] for epoch in range(10, 40)], c='blue', label='Basic')

plt.legend(bbox_to_anchor=(1.1, 0.5))
plt.show()
plt.figure(figsize=(10,5))
plt.plot([float(x) for x in hydra_metrics['val_vae_loss_1__amino_acc'].tolist()[10:40]],
         c='red', label='HydrAMP', linestyle='--')
plt.plot([float(x) for x in pepcvae_metrics['val_vae_loss_1__amino_acc'].tolist()[10:40]],
         c='orange', label='PepCVAE', linestyle='--')
plt.plot([float(x) for x in basic_metrics['val_vae_loss_1__amino_acc'].tolist()[10:40]],
         c='blue', label='Basic', linestyle='--')



plt.title('How many petides were susceptible to (relative) improvement out of 2404 known AMPs? ')
plt.plot([pos_final_results[models[0]][epoch]['before relative improvement']/2404 for epoch in range(10, 40)], c='red')
plt.plot([pos_final_results[models[1]][epoch]['before relative improvement']/2404 for epoch in range(10, 40)], c='orange')
plt.plot([pos_final_results[models[2]][epoch]['before relative improvement']/2404 for epoch in range(10, 40)], c='blue')

plt.legend(bbox_to_anchor=(1.1, 0.5))
plt.show()
plt.figure(figsize=(10,5))
plt.plot([float(x) for x in hydra_metrics['val_vae_loss_1__amino_acc'].tolist()[10:40]],
         c='red', label='HydrAMP', linestyle='--')
plt.plot([float(x) for x in pepcvae_metrics['val_vae_loss_1__amino_acc'].tolist()[10:40]],
         c='orange', label='PepCVAE', linestyle='--')
plt.plot([float(x) for x in basic_metrics['val_vae_loss_1__amino_acc'].tolist()[10:40]],
         c='blue', label='Basic', linestyle='--')


plt.title('How many petides were susceptible to (absolute) improvement out of 2404 known AMPs? ')
plt.plot([pos_final_results[models[0]][epoch]['before absolute improvement']/2404 for epoch in range(10, 40)], c='red')
plt.plot([pos_final_results[models[1]][epoch]['before absolute improvement']/2404 for epoch in range(10, 40)], c='orange')
plt.plot([pos_final_results[models[2]][epoch]['before absolute improvement']/2404 for epoch in range(10, 40)], c='blue')

plt.legend(bbox_to_anchor=(1.1, 0.5))
plt.show()<jupyter_output><empty_output><jupyter_text>Model selection<jupyter_code>def choose_best_epoch(model):
    model_metrics = pd.read_csv(f'../models/{model}/metrics.csv')
    good_epochs = model_metrics.iloc[10:40][model_metrics['val_vae_loss_1__amino_acc'].astype(float)  > 0.95].epoch_no.tolist()
    improved_peptides = [pos_final_results[model][epoch]['before relative improvement']/2404 + \
                         neg_final_results[model][epoch]['before absolute improvement']/2223 \
                         for epoch in good_epochs]
    return good_epochs[np.argmax(improved_peptides)], np.max(improved_peptides)  
best_epochs = {model: [] for model in models}

for model in models:
    best_epochs[model] = choose_best_epoch(model)

best_epochs
ax = sns.barplot(
    x=model_labels,
    y=[
        pos_final_results[model][int(best_epochs[model][0])]['before relative improvement']/2404 + \
        neg_final_results[model][int(best_epochs[model][0])]['before absolute improvement']/2223 \
        for model in models
    ]
)

ax.set_title('VALIDATION SET\n % of relatively improved positives + % of absolutely improved negatives')
ax.set_xticklabels(model_labels, rotation=90)
plt.show()
    


metrics_to_consider = [
    'before relative improvement',
    'before absolute improvement',
    'relative improvement',
    'absolute improvement',
]

metrics_labels = [
    'How many petides were susceptible to (relative) improvement?',
    'How many petides were susceptible to (absolute) improvement?',
    'Number of uniquely generated peptides during relative improvement procedure (64 attempts per peptide)',
    'Number of uniquely generated peptides during absolute improvement procedure (64 attempts per peptide)',
]
for metric, metric_label in zip(metrics_to_consider, metrics_labels):    
    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 4), sharex=True)
    plt.suptitle(metric_label, y=1.1)
    sns.barplot(x=model_labels, y=[pos_final_results[model][int(best_epochs[model][0])][metric] for model in models], ax=ax1)
    sns.barplot(x=model_labels, y=[neg_final_results[model][int(best_epochs[model][0])][metric] for model in models], ax=ax2)

    ax1.set_title('2404 positives (validation set)')
    ax2.set_title('2223 negatives (validation set)')
    ax1.set_xticklabels(model_labels, rotation=90)
    ax2.set_xticklabels(model_labels, rotation=90)
    plt.show()<jupyter_output><empty_output><jupyter_text>Test set<jupyter_code>best_epochs = {
    'HydrAMP': 37,
    'PepCVAE': 35,
    'Basic': 15,
}
pos = np.vstack([amp_x_test[amp_y_test == 1], mic_x_test[mic_y_test < 1.5]])
neg = np.vstack([amp_x_test[amp_y_test == 0], mic_x_test[mic_y_test > 1.5]])

print(pos.shape, neg.shape)

pos_amp = amp_classifier_model.predict(pos, verbose=1).reshape(len(pos))
neg_mic = mic_classifier_model.predict(neg, verbose=1).reshape(len(neg))
neg_amp = amp_classifier_model.predict(neg, verbose=1).reshape(len(neg))
pos_mic = amp_classifier_model.predict(pos, verbose=1).reshape(len(pos))

pos = np.vstack([pos] * 64).reshape(-1, 25)
pos_amp = np.vstack([pos_amp] * 64).reshape(-1, 1)
pos_mic = np.vstack([pos_mic] * 64).reshape(-1, 1)
neg = np.vstack([neg] * 64).reshape(-1, 25)
neg_amp = np.vstack([neg_amp] * 64).reshape(-1, 1)
neg_mic = np.vstack([neg_mic] * 64).reshape(-1, 1)
final_pos_results = {}
final_neg_results = {}
for model in tqdm(models):
    epoch = int(best_epochs[model])
    AMPMaster = bms.load_model(f'../models/{model}/{epoch}')
    encoder_model =  AMPMaster.encoder(input_to_encoder)
    decoder_model = AMPMaster.decoder(input_to_decoder)
    if model in ['PepCVAE', 'Basic']:
        new_act = layers.TimeDistributed(
            layers.Activation(activations.softmax),
            name='decoder_time_distribute_activation')
        decoder_model.layers.pop()
        x = new_act(decoder_model.layers[-1].output)
        decoder_model = Model(input=decoder_model.input, output=[x]) 
    final_pos_results[model] = improve(pos, model, epoch, 'pos')
    final_neg_results[model] = improve(neg, model, epoch, 'neg')
    
dump(final_pos_results, f'../results/improvement_PosTest.joblib') 
dump(final_neg_results, f'../results/improvement_NegTest.joblib') 


pos_final_results = {models: {} for models in models}
neg_final_results = {models: {} for models in models}
for model in models:
        pos_final_results[model]['relative improvement'] = np.unique(
            final_pos_results[model]['rel_improved'], axis=0).shape[0]
        pos_final_results[model]['absolute improvement'] = np.unique(
             final_pos_results[model]['abs_improved'], axis=0).shape[0]
        pos_final_results[model]['before relative improvement'] = np.unique(
            final_pos_results[model]['before_rel_improve'], axis=0).shape[0]
        pos_final_results[model]['before absolute improvement'] = np.unique(
            final_pos_results[model]['before_abs_improve'], axis=0).shape[0]
        
        neg_final_results[model]['relative improvement'] = np.unique(
            final_neg_results[model]['rel_improved'], axis=0).shape[0]
        neg_final_results[model]['absolute improvement'] = np.unique(
             final_neg_results[model]['abs_improved'], axis=0).shape[0]
        neg_final_results[model]['before relative improvement'] = np.unique(
            final_neg_results[model]['before_rel_improve'], axis=0).shape[0]
        neg_final_results[model]['before absolute improvement'] = np.unique(
            final_neg_results[model]['before_abs_improve'], axis=0).shape[0]   
ax = sns.barplot(
    x=model_labels,
    y=[
        pos_final_results[model]['before relative improvement']/1319 + \
        neg_final_results[model]['before absolute improvement']/1253 \
        for model in models])

ax.set_title('Before relative improvement (positives) + before absolute improvement (negatives)')
ax.set_xticklabels(model_labels, rotation=90)
plt.show()
    
for metric, metric_label in zip(metrics_to_consider, metrics_labels):    
    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 4), sharex=True)
    plt.suptitle(metric_label, y=1.1)
    sns.barplot(x=model_labels, y=[pos_final_results[model][metric] for model in models], ax=ax1)
    sns.barplot(x=model_labels, y=[neg_final_results[model][metric] for model in models], ax=ax2)

    ax1.set_title('1319 positives (test set)')
    ax2.set_title('1253 negatives (test set)')
    ax1.set_xticklabels(model_labels, rotation=90)
    ax2.set_xticklabels(model_labels, rotation=90)
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14863251573692207, 'fraction_numerical': 0.026372910787931407, 'mean_word_length': 4.399648403164371, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '10806', 'license': 'MIT', 'n_tokens_mistral': 6778, 'n_tokens_neox': 6353, 'n_words': 1077, 'path': 'scripts/model_selection.ipynb', 'repo_name': 'szczurek-lab/hydramp', 'size': '18428'}