{"cells":[{"metadata":{},"cell_type":"markdown","source":"**Transfer learning framework with keras starter kernel by alex kharin**\nHere is the starting kernel used Xception, darknet or any other ConvNN for prediction of the cell treatment.\nHere is pipeline, where you can play with different base models for feature extraction and different classifier models architectures. Feature extraction and further classifier NN are different networks, which can imporove the development speed. Do not hesitate to change Createmodel() function and basemodel. Other thing to play with is dimensionality reduction ways.\n\n\nV34 Changed resolution of the input images to original, changing pretrained backbone input shape, separate models for each celltype, use of preextracted features ()\n\nto add: kfolds\n\nV33\nminor changes in top layer architecture; addition of LR exponential decay\n    \nversion 31 - Fixed mistake in path construction for control samples\n\nversion 30 - Features to dense NN are now the differences in features extracted by base_model in train images and corresponding negative controls\ngeneral pipeline of the kernel\nV30 pipeline:\n    1. Create dataframe with all  train samples and filepaths\n    2. Split on train/val\n    3. Make linear dimensionality reduction (each channel is linear combination of initial ones, weights can be obtained by PCA or linear autoencoder)\n    4. Create generator of extracted features from images using pretrained models (Xception, darkent or smth else) and quantile extraction\n    5. Extract features and (optionally) save them for faster development\n    6. (since version 30) Calculate difference between features extracted from train images and corresponding negative control images\n    7. Create and train simple dense nn model with input shape corresponding to output of feature extractor\n    8. Save model (can be done using checkpoints) \n    9. Predict on test data with taking in account test-negative controls\n\nversion 29 - First satisfactory version according to general pipeline\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"#!pip install keras-mxnet kerascv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#import  kerascv\n#from kerascv.model_provider import get_model as kecv_get_model","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"#import gcimport gc\nimport os\nimport pickle\nimport random\nimport time\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport PIL\nfrom PIL import Image, ImageOps\nimport sklearn\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.decomposition import PCA\nimport keras\nimport matplotlib\nimport matplotlib.pyplot as plt\nfrom keras.models import Model\nimport csv\nfrom keras.layers import GlobalAveragePooling2D, GlobalMaxPooling2D, Dense\nimport tensorflow as tf\nimport time\nfrom keras.engine.topology import Layer\nimport keras.backend as K","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def constructpath(row,train=True):\n    '''\n    function to construct filepath for a single specific row\n    '''\n    ans=[]\n    if train==True:\n        path='../input/train/'\n    else:\n        path='../input/test/'\n    experiment, plate, well =row['id_code'].split('_')\n    path=path+str(experiment+'/'+'Plate'+plate+'/'+ well)\n    for site in ['_s1','_s2']:\n        for channel in ['_w1','_w2','_w3','_w4','_w5','_w6']:\n            fp=path+site+channel+'.png'\n            #if os.path.exists(fp):\n            ans.append(fp)\n    ans.sort()\n    ans2=ans[6:]\n    ans=ans[:6]\n    return (tuple(ans),tuple(ans2))\n\ndef reduce_dims(batch, weights):\n    return np.dot(batch,weights[:,:])\n\ndef picgenerator3D(\n    trainDS2,Weights, basemodel,graph, site,batchsize=8, \n                 channelnumlist=[0,1,2,3,4,5],\n                 shuffle=True,autoenc=False, \n                 test=False,\n    ):\n    trainDS=trainDS2.copy().reset_index(drop=True)\n    shape=512\n    try: shape=int(basemodel.input.get_shape()[1])\n    except:shape=512\n    while True:\n        if shuffle==True:\n            trainDS = trainDS.sample(frac=1).reset_index(drop=True)\n        batch=[]\n        ans=np.zeros((batchsize)).astype(int)\n        for i, row in trainDS.iterrows():\n            if test==False:\n                ans[i%batchsize]=row['sirna']\n            imageC=[]\n            for channel in channelnumlist:\n                imageC.append(Image.open(row['pathes'][site][channel]).resize((shape,shape)))\n            batch.append(np.stack(imageC,axis=-1))\n            if (i+1)%batchsize==0:\n                with graph.as_default():\n                    out=basemodel.predict(reduce_dims(np.stack(batch)/255, Weights))\n                    out2=(np.quantile(out,0.9,axis=(1,2))+np.quantile(out,0.95,axis=(1,2)))/2\n                if autoenc==True:\n                    yield out2,out2\n                else:\n                    if test==False:\n                        yield out2, ans\n                    if test==True:\n                        yield out2\n                batch=[]\n                ans=np.zeros((batchsize))\n        if autoenc==True:\n            with graph.as_default():\n                out=basemodel.predict(reduce_dims(np.stack(batch)/255, Weights))\n                out2=(np.quantile(out,0.9,axis=(1,2))+np.quantile(out,0.95,axis=(1,2)))/2\n            yield out2,out2\n        else:\n            with graph.as_default():\n                out=basemodel.predict(reduce_dims(np.stack(batch)/255, Weights))\n                out2=(np.quantile(out,0.9,axis=(1,2))+np.quantile(out,0.95,axis=(1,2)))/2\n            if test==False:\n                yield out2, ans[:len(batch)]\n            if test==True:\n                yield out2\n                \ndef Createmodel(basemodel, layers,reg=0.0,lr=0.001,regllast=0.0,arcface=False):\n    shape=int(basemodel.output.shape[-1])\n    inp=keras.layers.Input(shape=(shape,))\n    label=keras.layers.Input(shape=(1108,))\n    if layers[0]!=0:\n        x=keras.layers.Dense(layers[0],activation='elu',kernel_regularizer=keras.regularizers.l2(reg),bias_regularizer=keras.regularizers.l2(reg))(inp)\n        x=keras.layers.Dropout(0.2)(x)\n        for obj in layers[1:]:\n            x=keras.layers.Dense(obj,activation='relu',kernel_regularizer=keras.regularizers.l2(reg),bias_regularizer=keras.regularizers.l2(reg))(x)\n            x=keras.layers.Dropout(0.2)(x)\n    #x=keras.layers.Dense(64,activation='relu',kernel_regularizer=keras.regularizers.l2(0.0000),bias_regularizer=keras.regularizers.l2(0.00000))(x)\n    #x=keras.layers.Dropout(01)(x)\n    if layers[0]==0:\n        x=keras.layers.Dropout(0.2)(inp)\n    if arcface==True:out=ArcFace(1108, 30, 0.5)([x,label]) \n    if arcface==False:\n        x=keras.layers.BatchNormalization()(x)\n        x=keras.layers.Dense(1108,activation='elu',kernel_regularizer=keras.regularizers.l2(regllast),bias_regularizer=keras.regularizers.l2(regllast))(x)\n        out=keras.layers.Activation('softmax')(x)\n    if arcface==True:model=Model(input=[inp,label],outputs=out)\n    if arcface==False:model=Model(input=inp,outputs=out)\n    model.compile(optimizer=keras.optimizers.Adam(lr=lr), loss='categorical_crossentropy',metrics=['categorical_accuracy'])\n    return model\n\ndef featurecollector(DF, basemodel,weights, batchsize,site=1,test=False):\n    inittime=time.time()\n    generator=picgenerator3D(\n    DF,weights,basemodel,graph,site,\n    batchsize=BATCHSIZE,channelnumlist=[0,1,2,3,4,5], \n    autoenc=False,\n    test=test,\n    shuffle=False\n        )\n    toret=np.zeros((len(DF),int(basemodel.output.shape[-1])))\n    ysez=np.zeros(len(DF))\n    #print (toret.shape[0])\n    for it in range(int(np.ceil(len(DF)/batchsize))):\n        A=next(generator)\n        if test==False:\n            toret[it*batchsize:it*batchsize+A[0].shape[0],:]=A[0]\n            ysez[it*batchsize:it*batchsize+A[1].shape[0]]=A[1]\n            timefrombeg=time.time()-inittime\n        else:\n            toret[it*batchsize:it*batchsize+A.shape[0],:]=A\n            #ysez[it*batchsize:it*batchsize+A[1].shape[0]]=A[1]\n            timefrombeg=time.time()-inittime\n        if it%30==0 and it>0: print ('remains', timefrombeg/it*(int(np.ceil(len(DF)/batchsize)-it)), ' s')\n    return toret, ysez\ndef schedule(epoch,lr):\n    return lr*0.9","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras import backend as K\nfrom keras.layers import Layer\nfrom keras import regularizers\n\nimport tensorflow as tf\nclass ArcFace(Layer):\n    def __init__(self, n_classes=10, s=30.0, m=0.50, regularizer=None, **kwargs):\n        super(ArcFace, self).__init__(**kwargs)\n        self.n_classes = n_classes\n        self.s = s\n        self.m = m\n        self.regularizer = regularizers.get(regularizer)\n\n    def build(self, input_shape):\n        super(ArcFace, self).build(input_shape[0])\n        self.W = self.add_weight(name='W',\n                                shape=(input_shape[0][-1], self.n_classes),\n                                initializer='glorot_uniform',\n                                trainable=True,\n                                regularizer=self.regularizer)\n\n    def call(self, inputs):\n        x, y = inputs\n        c = K.shape(x)[-1]\n        # normalize feature\n        x = tf.nn.l2_normalize(x, axis=1)\n        # normalize weights\n        W = tf.nn.l2_normalize(self.W, axis=0)\n        # dot product\n        logits = x @ W\n        # add margin\n        # clip logits to prevent zero division when backward\n        theta = tf.acos(K.clip(logits, -1.0 + K.epsilon(), 1.0 - K.epsilon()))\n        target_logits = tf.cos(theta + self.m)\n        # sin = tf.sqrt(1 - logits**2)\n        # cos_m = tf.cos(logits)\n        # sin_m = tf.sin(logits)\n        # target_logits = logits * cos_m - sin * sin_m\n        #\n        logits = logits * (1 - y) + target_logits * y\n        # feature re-scale\n        logits *= self.s\n        out = tf.nn.softmax(logits)\n\n        return out\n\n    def compute_output_shape(self, input_shape):\n        return (None, self.n_classes)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weights=np.array([[ 0.36013198, -0.        ,  0.10044181],\n       [ 0.04300603,  0.27274185,  0.3570403 ],\n       [ 0.14618097,  0.20497076, -0.        ],\n       [ 0.19643596,  0.14860706,  0.50436   ],\n       [ 0.18471946,  0.22273754, -0.        ],\n       [ 0.0695256 ,  0.15094277,  0.03815791]], dtype=float)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"HERE YOU CAN CHANGE TO OTHER BASE MODEL FROM https://pypi.org/project/kerascv/ or from keras application and choose output layer of shape (None,a,b,feature_number). "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Data reading and preprocessing\nBATCHSIZE=32\nDF=pd.read_csv('../input/recursion-cellular-image-classification/train.csv')\nDF['pathes']=DF.apply(lambda x: constructpath(x),axis=1)\nDF['celltype']=DF['experiment'].apply(lambda x: x.split('-')[0])\nDF_train_control=pd.read_csv('../input/recursion-cellular-image-classification/train_controls.csv')\nDF_train_control['pathes']=DF_train_control.apply(lambda x: constructpath(x),axis=1)\nDF_test_control=pd.read_csv('../input/recursion-cellular-image-classification/test_controls.csv')\nDF_test_control['pathes'] = DF_test_control.apply(lambda x: constructpath(x, train = False), axis = 1)\n\n#get the samples by celltype\nidxHEPG2=list(DF[DF['celltype']=='HEPG2'].index)\nidxHUVEC=list(DF[DF['celltype']=='HUVEC'].index)\nidxRPE=list(DF[DF['celltype']=='RPE'].index)\nidxU2OS=list(DF[DF['celltype']=='U2OS'].index)\n\nsubmittionDF=pd.read_csv('../input/recursion-cellular-image-classification/test.csv')\nsubmittionDF['pathes']=submittionDF.apply(lambda x: constructpath(x,train=False),axis=1)\nsubmittionDF['celltype']=submittionDF['experiment'].apply(lambda x: x.split('-')[0])\n\nidxHEPG2test=submittionDF[submittionDF['celltype']=='HEPG2'].index\nidxHUVECtest=submittionDF[submittionDF['celltype']=='HUVEC'].index\nidxRPEtest=submittionDF[submittionDF['celltype']=='RPE'].index\nidxU2OStest=submittionDF[submittionDF['celltype']=='U2OS'].index\n\n\n#change INPUT SHAPE from (224,224) to (512,512)\n\nbasemodel=keras.applications.xception.Xception(include_top=False, weights='imagenet', input_tensor=None, input_shape=None, pooling=None, classes=1108)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Extract features using backbone and save them\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n'''\nglobal graph\ngraph = tf.get_default_graph() \n\n#EXTRACT FEATURES from pictures using model it outputs some logs with time remain to the feature extraction finish\n\ntraincontrols1=featurecollector(DF_train_control,basemodel, weights,BATCHSIZE,site=0)\ntraincontrols2=featurecollector(DF_train_control,basemodel, weights,BATCHSIZE,site=1)\ntraincontrols=(np.mean([traincontrols1[0],traincontrols2[0]],axis=0), traincontrols1[1])\nMODELINPS1=featurecollector(DF,basemodel, weights,BATCHSIZE,site=0)\n#pickle.dump(MODELINPS1, open(\"firstsitesfeatures\",\"wb\"))\nMODELINPS2=featurecollector(DF,basemodel, weights,BATCHSIZE,site=1)\n#pickle.dump(MODELINPS2, open(\"secondsitesfeatures\",\"wb\"))\n# calculate difference in features compared to control\nMODELINPS_1=np.zeros((MODELINPS1[0].shape))\nMODELINPS_2=np.zeros((MODELINPS2[0].shape))\nfor i,row in DF.iterrows():\n    ind = DF_train_control[DF_train_control['experiment']==row['experiment']][DF_train_control['plate']==row['plate']][DF_train_control['well_type']=='negative_control'].index[0]\n    MODELINPS_1[i,:]=MODELINPS1[0][i,:]-traincontrols[0][ind,:]\n    MODELINPS_2[i,:]=MODELINPS2[0][i,:]-traincontrols[0][ind,:]\n    \n#reading and preprocessing the test data, creating test features on site 0 and site 1 \n\n\nfeaturestest1=featurecollector(submittionDF,basemodel,weights,BATCHSIZE,site=0,test=True)\nfeaturestest2=featurecollector(submittionDF,basemodel,weights,BATCHSIZE,site=1,test=True)\ntestcontrols1=featurecollector(DF_test_control,basemodel, weights,BATCHSIZE,site=0)\ntestcontrols2=featurecollector(DF_test_control,basemodel, weights,BATCHSIZE,site=1)\nfeaturestest_1=np.zeros((featurestest1[0].shape))\nfeaturestest_2=np.zeros((featurestest2[0].shape))\ntestcontrols=(np.mean([testcontrols1[0],testcontrols2[0]],axis=0),testcontrols1[1])\nfor i,row in submittionDF.iterrows():\n    ind = DF_test_control[DF_test_control['experiment']==row['experiment']][DF_test_control['plate']==row['plate']][DF_test_control['well_type']=='negative_control'].index[0]\n    featurestest_1[i,:]=featurestest1[0][i,:]-testcontrols[0][ind,:]\n    featurestest_2[i,:]=featurestest2[0][i,:]-testcontrols[0][ind,:]\n    \nXsez=np.concatenate([MODELINPS_1,MODELINPS_2],axis=0)\nYsez1=np.concatenate([MODELINPS1[1],MODELINPS2[1]],axis=0)\npickle.dump(Xsez, open(\"Xsez\",\"wb\"))\npickle.dump(Ysez1, open(\"Ysez1\",\"wb\"))\n\npickle.dump(featurestest_1, open('Xtest1','wb'))\npickle.dump(featurestest_2, open('Xtest2','wb'))\n\n'''\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"os.listdir('../input/xception-cells-features-extracted/')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xsez=pickle.load(open('../input/xception-cells-features-extracted/Xsez','rb'))\nYsez1=pickle.load(open('../input/xception-cells-features-extracted/Ysez1','rb'))\nXtest1=pickle.load(open('../input/xception-cells-features-extracted/Xtest1','rb'))\nXtest2=pickle.load(open('../input/xception-cells-features-extracted/Xtest2','rb'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Ysez=keras.utils.to_categorical(Ysez1, num_classes=1108)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Here and in Createmodel() function you can play with classification model\nmodels={}\nfor modelname in ['modelHEPG2','modelHUVEC','modelRPE','modelU2OS']:\n    models[modelname]=[\n    Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False),\n    Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False),\n    Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False),\n    Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False),\n    ]\n'''modelHEPG2=Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False)\nmodelHUVEC=Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False)\nmodelRPE=Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False)\nmodelU2OS=Createmodel(basemodel,layers=[128,128],reg=0.01,lr=0.001,regllast=0.001,arcface=False)'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#models['modelHEPG2'][0].summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(123)\nnp.random.shuffle(idxHEPG2)\nnp.random.shuffle(idxHUVEC)\nnp.random.shuffle(idxRPE)\nnp.random.shuffle(idxU2OS)\nindicies={\n    'modelHEPG2':idxHEPG2,\n    'modelHUVEC':idxHUVEC,\n    'modelRPE':idxRPE,\n    'modelU2OS':idxU2OS\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kf=sklearn.model_selection.KFold(n_splits=4)\nhist={}\nfor modelname in ['modelHEPG2','modelHUVEC','modelRPE','modelU2OS']: \n    hist[modelname]=[0,1,2,3]\nfor modelname in ['modelHEPG2','modelHUVEC','modelRPE','modelU2OS']:\n    i=0\n    for train_index, val_index in kf.split(indicies[modelname]):\n        hist[modelname][i]=models[modelname][i].fit(\n            Xsez[train_index],Ysez[train_index],\n            epochs=50, \n            batch_size=BATCHSIZE,\n            validation_data=(Xsez[val_index],Ysez[val_index]),\n            callbacks=[keras.callbacks.ModelCheckpoint(modelname+str(i)+'.hdf5', monitor='val_categorical_accuracy', verbose=0, save_best_only=True),keras.callbacks.LearningRateScheduler(schedule)]\n        )\n        i=i+1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for modelname in ['modelHEPG2','modelHUVEC','modelRPE','modelU2OS']:\n    for i in range(4):\n        models[modelname][i]=keras.models.load_model(modelname+str(i)+'.hdf5')\n#np.max(histHEPG2.history['val_categorical_accuracy'])\n#np.max(histHEPG2.history['categorical_accuracy'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#make predictions on site 0 and 1\nanswer1=np.zeros((len(Xtest1),1108))\nanswer2=np.zeros((len(Xtest2),1108))\nfor ind, modelname in zip([idxHEPG2test,idxHUVECtest,idxRPEtest,idxU2OStest],['modelHEPG2','modelHUVEC','modelRPE','modelU2OS']):\n    for i in range(4):\n        answer1[ind]=answer1[ind]+models[modelname][i].predict(Xtest1[ind],verbose=1)\n        answer2[ind]=answer2[ind]+models[modelname][i].predict(Xtest2[ind],verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create submission file\nsubmittionDF['sirna']=np.argmax(answer1+answer2, axis=1)\nsubmittionDF.to_csv('submission.csv', sep=',', columns=['id_code', 'sirna'], header=True, index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}