{"cells":[{"metadata":{},"cell_type":"markdown","source":"Hey! I'm interested in joining a team. I'm currently working as a programmer at the Howard Hughes Medical Institute at Janelia, processing large-scale neural calcium imaging recordings https://www.janelia.org/lab/pachitariu-lab I have a master's degree in Applied Mathematics and I've worked as a data scientist in a team focused on big data. I've also done an internship in applying convolutional neural networks on satellite images. My github is here https://github.com/mariakesa If you would be interested in taking me into your team shoot me an email at maria.kesa@gmail.com\n\nWavelets are signal transform methods that can find a sparse representation of a signal based on it's localized properties. The Discrete Wavelet Transform is a way to analyze a signal at multiple scales (multiple frequency bands) The outputs giving the detail coefficients (from the high-pass filter) and approximation coefficients (from the low-pass) (Wikipedia). For one transform operation there are three high pass coefficient sets and one low pass coefficient set. Here's an excellent article on using wavelets in machine learning (for example with convolutional neural networks  http://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/\n\nHere I sample a random image and make a widget to plot all of the wavelet bases approximations to the image implemented in the pytorch_wavelets library. \n\nThe plan is to combine this notebook with the previous notebook on PCA https://www.kaggle.com/mariakesa/pcaforvisualizingbatcheffects by computing PCA on wavelet coeffients extracted images from different batches to see if the confounding variables appear at different scales, e.g. would the batch effects be more clearly separable in high pass or low pass regimes (the plots would be analogous to the plots in the notebook pointed out earlier, but this time we would do PCA on the wavelet coefficients).\n\nFinally, it would be interesting to try to implement a combination of convolutional neural networks and wavelets, for example see https://arxiv.org/abs/1805.08620"},{"metadata":{},"cell_type":"markdown","source":"To enable ipywidgets in a jupyter notebook you have to run 'jupyter nbextension enable --py widgetsnbextension' in the terminal.Pip install ipywidgets via 'pip install ipywidgets'.\n\nThe pytorch_wavelets library has to be cloned from github for pip installation. The visualizations with the widgets are pretty cool. It might be worth it to download the notebook and run it locally:-)\n\nNB! Widgets can be moved one sample at a time by pressing the left and right arrow keys after clicking on the scroll bar!"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Instructions for installing pytorch_wavelets are here https://pytorch-wavelets.readthedocs.io/en/latest/readme.html\n#!pip install ipywidgets\n\n#!pip install PyWavelets\nimport ipywidgets as widgets\nimport numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nfrom sklearn.decomposition import PCA\n\nfrom PIL import Image\n\nimport torch\nimport torch.nn as nn\nimport torch.utils.data as D\nimport torch.nn.functional as F\n\nimport torchvision\nfrom torchvision import transforms as T\n\nfrom tqdm import tqdm\n\nimport matplotlib.colors as colors\nfrom matplotlib import cm\n\nfrom pytorch_wavelets import DWTForward\nimport pywt\nfrom torch.utils.data import RandomSampler","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path=\"../input\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#ImagesDS was taken from https://www.kaggle.com/leighplt/densenet121-pytorch\nclass ImagesDS(D.Dataset):\n    def __init__(self, csv_file, img_dir, mode='train', site=1, channels=[1,2,3,4,5,6]):\n        \n        df = pd.read_csv(csv_file)\n        self.records = df.to_records(index=False)\n        self.channels = channels\n        self.site = site\n        self.mode = mode\n        self.img_dir = img_dir\n        self.len = df.shape[0]\n        \n    @staticmethod\n    def _load_img_as_tensor(file_name):\n        with Image.open(file_name) as img:\n            return T.ToTensor()(img)\n\n    def _get_img_path(self, index, channel):\n        experiment, well, plate = self.records[index].experiment, self.records[index].well, self.records[index].plate\n        return '/'.join([self.img_dir,self.mode,experiment,f'Plate{plate}',f'{well}_s{self.site}_w{channel}.png'])\n        \n    def __getitem__(self, index):\n        paths = [self._get_img_path(index, ch) for ch in self.channels]\n        img = torch.cat([self._load_img_as_tensor(img_path) for img_path in paths])\n        \n        if self.mode == 'train':\n            return img, self.records[index].sirna\n        else:\n            return img, self.records[index].id_code\n\n    def __len__(self):\n        \"\"\"\n        Total number of samples in the dataset\n        \"\"\"\n        return self.len\n\nclass ProcessData():\n    def __init__(self, data_path):\n        self.data_path=data_path\n        self.batch_names={'HEPG2':['HEPG2-01','HEPG2-02','HEPG2-03','HEPG2-04','HEPG2-05','HEPG2-06','HEPG2-07'],\n                          'HUVEC':['HUVEC-01','HUVEC-02','HUVEC-03','HUVEC-04','HUVEC-05','HUVEC-06','HUVEC-07',\n                                  'HUVEC-08','HUVEC-09','HUVEC-10','HUVEC-11','HUVEC-12','HUVEC-13','HUVEC-14',\n                                  'HUVEC-15','HUVEC-16'],\n                          'RPE':['RPE-01','RPE-02','RPE-03','RPE-04','RPE-05','RPE-06','RPE-07'],\n                          'U2OS':['U2OS-01','U2OS-02','U2OS-03']                            \n        }\n\n    def create_random_loader(self, nr_of_samples):\n        dataset= ImagesDS(data_path+'/train.csv', self.data_path)\n        sampler = RandomSampler(np.arange(nr_of_samples, dtype=np.int64))\n        loader = torch.utils.data.DataLoader(dataset, batch_size=1, sampler=sampler, num_workers=2)\n        return loader\n    \n    def create_loader_for_pca(self,nr_of_samples,cell_type):\n        df=pd.read_csv(self.data_path+'/train.csv')\n        if cell_type=='HEPG2':\n            batches=self.batch_names['HEPG2']\n        if cell_type=='HUVEC':\n            batches=self.batch_names['HUVEC']\n        if cell_type=='RPE':\n            batches=self.batch_names['RPE']\n        if cell_type=='U2OS':\n            batches=self.batch_names['U2OS']\n        df_lst=[]\n        for batch in batches:\n            ind=df['experiment']==batch\n            sub=np.array(df[ind].index)\n            sub_=list(df[ind].index)\n            nr_of_samples_in_batch=len(sub_)\n            generate_random_numbers=np.random.randint(0,nr_of_samples_in_batch,nr_of_samples)\n            sub=sub[generate_random_numbers]\n            df_lst=df_lst+list(sub)\n        df_=df.loc[df_lst,:]\n        df_.to_csv(data_path+'/train_pca.csv')\n        dataset = ImagesDS(data_path+'/train_pca.csv', self.data_path)\n        loader = D.DataLoader(dataset, batch_size=1, shuffle=False, num_workers=1)\n        return loader\n        \n    def flatten_data_for_PCA(self,cell_type,nr_of_samples):\n        loader=self.create_loader_for_pca(nr_of_samples,cell_type)\n        arr=torch.zeros((1572864,1))\n        for x, y in tqdm(loader):\n            #Flatten into 1D vector of features for PCA\n            x=x.flatten().view(1572864,1)\n            arr=torch.cat((arr,x),dim=1)\n        return arr[:,1:].numpy()\n    \n    def PCA_for_batches(self,cell_type,n_components,nr_of_samples):\n        '''\n        nr of samples denotes how many samples to take from each batch. \n        '''\n        arr=self.flatten_data_for_PCA(cell_type,nr_of_samples)\n        pca=PCA(n_components=n_components)\n        tr=pca.fit_transform(arr.T)\n        plt.plot(np.cumsum(pca.explained_variance_ratio_))\n        plt.title(cell_type)\n        plt.xlabel('number of components')\n        plt.ylabel('cumulative explained variance')\n        plt.show()\n        labels=[]\n        lab_dict={}\n        colors_=range(0,len(self.batch_names[cell_type]))\n        cmp = cm.get_cmap('viridis',len(self.batch_names[cell_type]))\n        for j in range(0,len(self.batch_names[cell_type])):\n            labels=labels+[j]*nr_of_samples\n            lab_dict[self.batch_names[cell_type][j]]=j\n        labels=np.array(labels)\n        for g in np.unique(self.batch_names[cell_type]):\n            ix=np.where(labels==lab_dict[g])\n            plt.scatter(tr[:,0].flatten()[ix], tr[:,1].flatten()[ix], cmap=cmp, label=g)\n        plt.legend()\n        plt.title(cell_type)\n        plt.xlabel('component 1')\n        plt.ylabel('component 2')\n        plt.show()\n        \n    def take_all_wavelet_transforms_of_a_single_image(self):\n        wavelets=pywt.wavelist()\n        print(wavelets)\n        self.wavelets_=[]\n        excluded_continuous=['cgau1','cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7','cgau8','cmor','fbsp','gaus1', \n                             'gaus2', 'gaus3', 'gaus4', 'gaus5', 'gaus6', 'gaus7', 'gaus8','mexh','morl','shan']\n        for wavelet in wavelets:\n            if wavelet not in excluded_continuous:\n                self.wavelets_.append(wavelet)\n        loader=self.create_random_loader(nr_of_samples=1)\n        self.Yl_lst=[]\n        self.Yh0_lst=[]\n        self.Yh1_lst=[]\n        self.Yh2_lst=[]\n        for x,y in loader:\n            for wavelet in self.wavelets_:\n                xfm = DWTForward(J=1, mode='symmetric', wave=wavelet)  # Accepts all wave types available to PyWavelets\n                Yl, Yh = xfm(x)\n                Yl=Yl.numpy()\n                Yl_=Yl[0,:,:,:]\n                self.Yl_lst.append(Yl_)\n                Yh=Yh[0]\n                Yh0=Yh[0,:,0,:,:].numpy()\n                Yh1=Yh[0,:,1,:,:].numpy()\n                Yh2=Yh[0,:,2,:,:].numpy()\n                self.Yh0_lst.append(Yh0)\n                self.Yh1_lst.append(Yh1)\n                self.Yh2_lst.append(Yh2)\n                \n    def plot_wavelet_Yl(self,x):\n        for j in range(0,6):\n            plt.imshow(self.Yl_lst[x][j,:,:])\n            plt.title(self.wavelets_[x]+', Channel '+str(j))\n            plt.show()\n    \n    def plot_wavelet_Yh0(self,x):\n        for j in range(0,6):\n            plt.imshow(self.Yh0_lst[x][j,:,:])\n            plt.title(self.wavelets_[x]+', Channel '+str(j))\n            plt.show()\n            \n    def plot_wavelet_Yh1(self,x):\n        for j in range(0,6):\n            plt.imshow(self.Yh1_lst[x][j,:,:])\n            plt.title(self.wavelets_[x]+', Channel '+str(j))\n            plt.show()\n        \n    def plot_wavelet_Yh2(self,x):\n        for j in range(0,6):\n            plt.imshow(self.Yh2_lst[x][j,:,:])\n            plt.title(self.wavelets_[x]+', Channel '+str(j))\n            plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"proc=ProcessData(data_path)\nproc.take_all_wavelet_transforms_of_a_single_image()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"widgets.interact(proc.plot_wavelet_Yl, x=(0,len(proc.wavelets_)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Check out channel 79, rbio3.1 Looks nice.\nwidgets.interact(proc.plot_wavelet_Yh0, x=(0,len(proc.wavelets_)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"widgets.interact(proc.plot_wavelet_Yh1, x=(0,len(proc.wavelets_)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"widgets.interact(proc.plot_wavelet_Yh2, x=(0,len(proc.wavelets_)))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.8"}},"nbformat":4,"nbformat_minor":1}