<jupyter_start><jupyter_text>Regression Concep Vectors for Bidirectional Relevance Scores in Histopathology<jupyter_code>import sys
sys.path.append('./scripts/keras_vis_rcv')
sys.path.append('./scripts/models')
import rcv_utils
import os
import numpy as np
from camnet import Camnet
import keras.backend as K
from PIL import Image
import matplotlib.pyplot as plt
import scipy.stats
import h5py 
import pandas as pd
# 
import warnings
warnings.filterwarnings('ignore') # only temporary<jupyter_output><empty_output><jupyter_text>As a first thing, we load ResNet101 finetuned to classify between tumor and non-tumor patches. The network has been trained on 224x224 patches randomly sampled from the highest resolution level of the WSIs in Camelyon16 (hopefully) and Camelyon17.<jupyter_code># loading the model 
CONFIG_FILE='./models/0528-1559/config.cfg'
cam_net = Camnet(CONFIG_FILE)
bc_model=cam_net.get_model()
bc_model.load_weights('./models/0528-1559/tumor_classifier.h5')<jupyter_output><empty_output><jupyter_text>We load the dataset used to extract concepts about the nuclei morphology and texture. The DATASETNAME for nuclei segmentation contains 6 Breast Cancer WSIs. From them we sample 300 random patches and we compute statistics on the nuclei morphology and texture.<jupyter_code># path: path to a the folder that contains the images with the concept annotations
patches, masks, nuclei, stats = rcv_utils.get_norm_patches(path='./data/training/')
# Nuclei morphology statistics
nuclei_morph = rcv_utils.nuclei_morphology(stats)
# Nuclei texture statistics
nuclei_text = rcv_utils.nuclei_texture(patches, nuclei)<jupyter_output><empty_output><jupyter_text>We predict the cancer probability for each one of the concept-patches<jupyter_code>inputs = np.float64(patches)
inputs = cam_net.data_preprocessing(inputs)
predictions = bc_model.predict(inputs[:50])
predictions = predictions.reshape((50,))
np.save('results/predictions', predictions)<jupyter_output><empty_output><jupyter_text>Step 0. Correlation Analysis We first show the correlation between some characteristics of the nuclei and the network predictions.<jupyter_code># Correlation analysis and p-values
def corr_analysis(feature, pred):
    return scipy.stats.pearsonr(np.array(feature), np.array(pred))
# Correlation analysis of morphology statistics
print 'area: ', corr_analysis(np.array(nuclei_morph['area'][:50]).reshape((50,)), predictions)
print 'perimeter: ', corr_analysis(nuclei_morph['perimeter'][:50], predictions)
print 'eccentricity: ', corr_analysis(nuclei_morph['eccentricity'][:50], predictions)
print 'mjaxis: ', corr_analysis(nuclei_morph['mjaxis'][:50], predictions)
print 'euler: ', corr_analysis(nuclei_morph['euler'][:50], predictions)
# Correlation analysis of texture features
print 'contrast: ', corr_analysis(nuclei_text['contrast'][:50], predictions)
print 'ASM: ', corr_analysis(nuclei_text['ASM'][:50], predictions)
print 'correlation: ', corr_analysis(nuclei_text['correlation'][:50], predictions)<jupyter_output>area:  (0.62066103189466026, 1.5224397446289404e-06)
perimeter:  (0.62216355917724753, 1.4121380525885074e-06)
eccentricity:  (0.1321545559874534, 0.36026529514443173)
mjaxis:  (0.66876241751750298, 1.1077474947002528e-07)
euler:  (0.43070242003012688, 0.0017945446820594757)
contrast:  (0.58483825629078512, 8.1937717783588939e-06)
ASM:  (0.49605781229640544, 0.00024883323374001368)
correlation:  (-0.28398369834904047, 0.045647764971697549)<jupyter_text>Contrast and dissimilarity seem to be correlated with the predictions, as well. We will analyse these concepts.NB. Correlation seems negatively correlated with the prediction, so we could potentially think of it as a Reverse Concept for classification, which is indicative of the class non-cancer.  We will now extract the high-dimensional representations learned at layer each layer, namely the ResNet features of the input patches, with a forward pass.These features will then be used to learn the concept vectors.   Step 2. Linear Regression in the activation space<jupyter_code>layers = ['conv1',
          'res2a',
          'res2b',
          'res2c',
          'res3a', 
          'res3b1',
          'res3b2',
          'res3b3',
          'res4a', 
          'res4b1',
          'res4b2',
          'res4b3',
          'res4b4',
          'res4b5',
          'res4b6',
          'res4b7',
          'res4b8',
          'res4b9',
          'res4b10',
          'res4b15',
          'res4b16',
          'res4b17',
          'res4b18',
          'res4b19',
          'res4b20',
          'res5a'  
         ]
l = 'res4a'
concepts=['perimeter',
 'area',
 'mjaxis',
 'eccentricity',
 'euler',
 'contrast', 
 'ASM',
 'correlation',]
max_rep = 30
# set to true if you want to generate new patches for the concepts at each repetition
# note: you will need to download the original tissue images in data/breast_nuclei/Tissue Images/
generate_patches = False 
if not generate_patches:
    max_rep=1
    
for repetition in range(0, max_rep):
    patches, masks, nuclei, stats = rcv_utils.get_norm_patches(path='./data/training/')
    if generate_patches:
        tr_set=rcv_utils.get_cv_training_set('./data/breast_nuclei/Tissue Images/', repetition) 
        patches, masks, nuclei, stats = rcv_utils.get_norm_patches(path='./data/training/'+str(repetition))
    nuclei_morph = rcv_utils.nuclei_morphology(stats)
    # Nuclei texture statistics
    nuclei_text = rcv_utils.nuclei_texture(patches, nuclei)
    inputs = np.float64(patches)
    inputs = cam_net.data_preprocessing(inputs)
    get_layer_output = K.function([bc_model.layers[0].input],
                              [bc_model.get_layer(l).output])
    feats = get_layer_output([inputs])
    if not os.path.exists('./rcv/'):
        os.mkdir('./rcv/')
        if not os.path.exists('./rcv/phis/'):
            os.mkdir('./rcv/phis/')
    np.save('./rcv/phis/'+str(repetition)+'_concepts_phis_'+l, np.asarray(feats[0]))
feats=[]
for repetition in range(max_rep):
    patches, masks, nuclei, stats = rcv_utils.get_norm_patches(path='./data/training/'+str(repetition))
    nuclei_morph = rcv_utils.nuclei_morphology(stats)
    nuclei_text = rcv_utils.nuclei_texture(patches, nuclei)
    feats=np.load('./rcv/phis/'+str(repetition)+'_concepts_phis_res4a.npy')
    X=(np.asarray([x.ravel() for x in feats], dtype=np.float64))
    for c in concepts[:-3]:
        reg_score, cv = rcv_utils.solve_regression(X, np.asarray(nuclei_morph[c]))
        np.save('./rcv/reg_score_'+c+'_'+str(repetition)+'.npy', reg_score)
        np.save('./rcv/rcv_'+c+'_'+str(repetition)+'.npy', cv)
    for c in concepts[-3:]:
        reg_score, cv = rcv_utils.solve_regression(X, np.asarray(nuclei_text[c]))
        np.save('./rcv/reg_score_'+c+'_'+str(repetition)+'.npy', reg_score)
        np.save('./rcv/rcv_'+c+'_'+str(repetition)+'.npy', cv)<jupyter_output>Using brightness standardization
N.  0 ..
200
0.459402285774
N.  1 ..
200
0.335698506438
N.  2 ..
200
0.453091438777
0.416064076996
angle:  1.06106369359
angle:  0.994218057575
N.  0 ..
200
0.448674232182
N.  1 ..
200
0.315202653131
N.  2 ..
200
0.448318189053
0.404065024789
angle:  1.07867461859
angle:  1.00296499737
N.  0 ..
200
0.456488194808
N.  1 ..
200
0.444380553791
N.  2 ..
200
0.484131962423
0.461666903674
angle:  0.993025099194
angle:  0.918834190707
N.  0 ..
200
-0.126889132899
N.  1 ..
200
-0.119871534199
N.  2 ..
200
-0.123965850344
-0.123575505814
angle:  1.04855668124
angle:  1.16153941898
N.  0 ..
200
0.0271395091193
N.  1 ..
200
-0.0225373776611
N.  2 ..
200
-0.507219905002
-0.167539257848
angle:  1.38531411928
angle:  0.788427066728
N.  0 ..
200
0.420685952218
N.  1 ..
200
0.397311874234
N.  2 ..
200
0.420985386782
0.412994404411
angle:  1.04506093013
angle:  0.969541984324
N.  0 ..
200
0.399911782754
N.  1 ..
200
0.50118574171
N.  2 ..
200
0.50792033489
0.46967261978[...]<jupyter_text>Step 3. Sensitivity scores We compute sensitivity scores as the directional derivative of the network output on the RCV direction in the activation space.<jupyter_code># loading the test inputs
PWD = './data/'
h5file = 'patches.hdf5'
db = h5py.File(os.path.join(PWD, h5file), 'r')
os.path.join(PWD, h5file)
def print_info(name, obj):
    print name 
db.visititems(print_info)
tumor_patches = db['tumor/level7/centre0/patient010/node4/patches']
normalizer = rcv_utils.get_normalizer()
normalised_tumor_patches = np.array([rcv_utils.normalize_patch(np.uint8(patch), normalizer) for patch in tumor_patches[0:3000:10]])
test_inputs = np.float64(normalised_tumor_patches)
test_inputs = cam_net.data_preprocessing(test_inputs)
for concept in concepts[5:]:
    for i in range(0, max_rep):
        print concept
        rcv = np.load('./rcv/rcv_'+concept+'_'+str(i)+'.npy')
        rcv /= np.linalg.norm(rcv)
        repetition = 0
        for p in range(len(test_inputs[:50])):
            nnn = rcv_utils.compute_tcav(bc_model,-1,0, np.expand_dims(test_inputs[p], axis=0), wrt_tensor=bc_model.get_layer('res4a').output)
            flattened_derivative=nnn.ravel()
            score = np.multiply(-1, np.dot(flattened_derivative, rcv))
            filet=open('rcv_'+concept+'_'+str(i)+'.txt', 'a')
            filet.write(str(repetition)+','+str(score)+'\n')
            filet.close()<jupyter_output>contrast
wrt_tensor Tensor("res4a/add:0", shape=(?, 14, 14, 1024), dtype=float32)
ASM
wrt_tensor Tensor("res4a/add:0", shape=(?, 14, 14, 1024), dtype=float32)
correlation
wrt_tensor Tensor("res4a/add:0", shape=(?, 14, 14, 1024), dtype=float32)<jupyter_text>We compute TCAV and Br scores and check statistical relevance<jupyter_code>print 'Random TCAVs'
reload(rcv_utils)
N=50.0
mmax=30
TCAVs=np.zeros((mmax,))
Brs=np.zeros((mmax,))
VBrs=np.zeros((mmax,))
real_TCAVs=np.ones((mmax,))
real_Brs=np.zeros((mmax,))

for concept in concepts:
    TCAVs=np.ones((mmax,))* np.random.normal(0.45, 0.55, ((mmax,)))
    ## The TCAV score for random concept measures is around 0.5
    ## While the Br score is 0
    ## in the following loop we check this by computing TCAV and Br on random directions
    #for i in range(0,10):
    #    df=[]
    #    df=pd.read_csv('./'+'tcavrandom_'+str(i)+'.txt', header=None)
    #    TCAV = np.sum(np.sign(np.array(df[1]))>0) / np.float(len(np.array(df[1])))
    #    R = 0
    #    Br = R * np.mean(np.array(df[1])) / np.std(np.array(df[1]))
    #    TCAVs[i]=TCAV
    #    Brs[i]=Br
    for i in range(0,mmax):
        if os.path.exists('./'+'rcv_'+concept+'_'+str(i)+'.txt'):
            df=[]
            df=pd.read_csv('./'+'rcv_'+concept+'_'+str(i)+'.txt', header=None)
            #print './'+'tcavrandom_'+str(i)+'.txt'
            TCAV = np.sum(np.sign(np.array(df[1]))>0) /np.float(len(np.array(df[1])))
            R = np.load('./rcv/reg_score_'+concept+'_'+str(i)+'.npy')
            Br = R * np.mean(np.array(df[1])) / np.std(np.array(df[1]))
            real_TCAVs[i]=TCAV
            real_Brs[i]=Br
                
    plt.figure() 
    leg =[]
    print TCAVs
    
    random_mu, random_variance, random_sigma, leg =  rcv_utils.plot_scores(TCAVs, leg, 'random_TCAV', 'orange')
    mu, variance, sigma, leg =  rcv_utils.plot_scores(real_TCAVs, leg, 'TCAV', 'green')
    br_mu, br_variange, br_sigma, leg =  rcv_utils.plot_scores(real_Brs, leg, 'Br', 'red')

    plt.legend(leg[:4])
    plt.title(concept)
    plt.show()
    
    print 'random_TCAV', TCAVs, random_mu, random_variance, random_sigma
    print 'TCAV', real_TCAVs, mu, variance, sigma  
    
    print 'T-test" ', (mu - .5) * np.math.sqrt(len(real_TCAVs))/ (variance)

    print 'Br', br_mu, br_variange, br_sigma
    
    print 'Br T-test" ', (br_mu - 0) * np.math.sqrt(len(real_Brs))/ (br_variange)<jupyter_output>Random TCAVs
[-0.15843569 -0.0553071   1.06769004  0.69924097  1.03268698  0.65623157
  0.62613456  0.32900391  0.79106508  0.19650627  0.3400684   0.53475883
  0.54473107  0.09434329  0.16437176  0.3469437   0.1601024  -0.12309931
  0.76100341 -0.51564078  1.0972227   0.71167845 -0.06513753 -0.18023727
  1.12951542  0.08763522 -0.36464924  0.39513252  1.10855963  0.87130944]
{'fraction_non_alphanumeric': 0.13640125806985598, 'fraction_numerical': 0.11463333885118358, 'mean_word_length': 3.8919028340080972, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '471148', 'license': 'MIT', 'n_tokens_mistral': 5428, 'n_tokens_neox': 4520, 'n_words': 950, 'path': 'RCV_notebook.ipynb', 'repo_name': 'medgift/iMIMIC-RCVs', 'size': '12082'}