{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Description\nThis kernel performs inference for [PANDA concat tile pooling starter](https://www.kaggle.com/iafoss/panda-concat-fast-ai-starter) kernel with use of multiple models and 8 fold TTA. Check it for more training details. The image preprocessing pipline is provided [here](https://www.kaggle.com/iafoss/panda-16x128x128-tiles).","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import cv2\nfrom tqdm import tqdm_notebook as tqdm\nimport fastai\nfrom fastai.vision import *\nimport os\n#from mish_activation import *\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport skimage.io\nimport numpy as np\nimport pandas as pd\n#sys.path.insert(0, '../input/semisupervised-imagenet-models/semi-supervised-ImageNet1K-models-master/')\n#from hubconf import *\nimport sys\nsys.path.append('/kaggle/input/panda-utils-cleaned/') #sys.path.append('/kaggle/input/panda-utils/')\nfrom model.model_config import ModelConfig\nfrom model.model_utils import load_models_from_dir","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA = '../input/prostate-cancer-grade-assessment/test_images'\nTEST = '../input/prostate-cancer-grade-assessment/test.csv'\nSAMPLE = '../input/prostate-cancer-grade-assessment/sample_submission.csv'\n\nbs = 2\nnworkers = 2\ncoef = [0.5, 1.5, 2.5, 3.5, 4.5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Change these settings","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model_dirs = ['../input/jj20200629', '../input/jj20200701/', '../input/jj20200703/', '../input/jj20200706']\nmodel_weights = np.array([3,4,2,1])\nuse_models = [[True,True,False,False], [True,True,False,False], [True,True,False,False], [True,False,False,False]]\ntile_size = [384,256,256,299]\nlevel=[1,1,1,1]\nmosaic_grid = [(4,6),(6,6),(6,6),(5,5)]\ntissue_th = (0.2,0.7)\nseed = 123\nsampling_method=['skeleton','skeleton','skeleton','skeleton']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The next cell will load your models and read additional settings from the model config.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# load config\nconfig = [ModelConfig.fromDir(model_dir) for model_dir in model_dirs]\n# load models\nmodels = [load_models_from_dir(model_dir, tile_list_input=False) for model_dir in model_dirs]\ncms = [[np.load(os.path.join(model_dir, f'cm_fold-{i}.npy')) for i in range(len(models[0]))] for model_dir in model_dirs]\n\nmodels = [[model for model, use_model in zip(_models,_use_models) if use_model] for _models,_use_models in zip(models,use_models)]\ncms = [[cm for cm, use_model in zip(_cms,_use_models) if use_model] for _cms,_use_models in zip(cms,use_models)]\n\nsz = [_config.getField('sz') for _config in config]\nN = [_config.getField('N') for _config in config]\nis_ordinal = [_config.getMetaField('is_ordinal')==True for _config in config]\nmodel_name = [_config.getField('model_name') for _config in config]\nregr = [\"regr\" in _model_name for _model_name in model_name]\nmean = [torch.tensor(np.array(_config.getField('mean')).astype(np.float32)) for _config in config]\nstd = [torch.tensor(np.array(_config.getField('std')).astype(np.float32)) for _config in config]\n#assert np.multiply(*mosaic_grid) == N, \"mosaic grid is different from the N read from the config\"\n\nfor i, _models in enumerate(models):\n    print(\"Loaded {0} {1} models\".format(len(_models), model_name[i]))\n    print(f'Mean {mean[i]} and std {std[i]}')\n    print(f'regression = {regr[i]}')\n    print(f'ordinal regression = {is_ordinal[i]}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(sz, tile_size)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## slide_preprocessing.py","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nfrom skimage.io import MultiImage\nfrom skimage.morphology import skeletonize\n#import maskslic as seg\n\nimport cv2\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nimport warnings\n\nVALID_SLIDE_EXTENSIONS = {'.tiff', '.mrmx', '.svs'}\n\n\n# ~~~~~~~~~~~~ Helper functions ~~~~~~~~~~~~\ndef generateMetaDF(data_dir, meta_fn:str='train.csv'):\n    '''\n        Makes a pandas.DataFrame of paths out of a directory including slides. Drop the `train.csv` in `data_dir`\n        and the script will also merge any meta data from there on `image_id` key.\n    '''\n    \n    \n    all_files = [path.resolve() for path in Path(data_dir).rglob(\"*.*\")]\n    slide_paths = [path for path in all_files if path.suffix in VALID_SLIDE_EXTENSIONS]\n    \n    if len(slide_paths)==0:\n        raise ValueError('No slides in `data_dir`=%s'%data_dir)\n    \n    data_df = pd.DataFrame({'slide_path':slide_paths})\n    data_df['image_id'] = data_df.slide_path.apply(lambda x: x.stem)\n    \n    slides = data_df[~data_df.image_id.str.contains(\"mask\")]\n    masks = data_df[data_df.image_id.str.contains(\"mask\")]\n    masks['image_id'] = masks.image_id.str.replace(\"_mask\", \"\")\n    masks.columns = ['mask_path', 'image_id']\n    \n    data_df = slides.merge(masks, on='image_id', how='left')\n    data_df['slide_path'] = data_df.slide_path.apply(lambda x: str(x) if not pd.isna(x) else None)\n    data_df['mask_path'] = data_df.mask_path.apply(lambda x: str(x) if not pd.isna(x) else None)\n    \n\n    ## Merge metadata\n    meta_csv = [file for file in all_files if file.name==meta_fn]\n    if meta_csv:\n        meta_df = pd.read_csv(str(meta_csv[0]))\n        data_df = data_df.merge(meta_df, on='image_id')\n    \n    return data_df\n\ndef tileClassification(tile, provider:str):\n    '''\n        Returns the cancer class of a tile based on majority vote of the tile's annotated pixels\n            0: background (non tissue) or unknown\n            1: benign tissue (stroma and epithelium combined)\n            2: cancerous tissue (stroma and epithelium combined)\n    '''\n\n    if provider == \"karolinska\":\n        '''\n        Karolinska: Regions are labelled. Valid values are:\n            0: background (non tissue) or unknown\n            1: benign tissue (stroma and epithelium combined)\n            2: cancerous tissue (stroma and epithelium combined)\n        '''\n        classes = {0:0, 1:1, 2:2}    \n        \n    elif provider == \"radboud\":\n        ''' \n        Radboud: Prostate glands are individually labelled. Valid values are:\n            0: background (non tissue) or unknown\n            1: stroma (connective tissue, non-epithelium tissue)\n            2: healthy (benign) epithelium\n            3: cancerous epithelium (Gleason 3)\n            4: cancerous epithelium (Gleason 4)\n            5: cancerous epithelium (Gleason 5)\n        '''\n        classes = {0:0, 1:1, 2:1, 3:2, 4:2, 5:2}\n    \n    tile = np.int32(tile)\n    counts = np.bincount(tile.reshape(-1,1)[:,0])\n\n    ## If only background, accept background (class=0) as the annotation\n    if len(counts)==1:\n        return 0\n\n    ## Otherwise take the second most common pixel as the annotation\n    max_annotation = np.argmax(counts[1:])       \n    \n    return classes[max_annotation+1]\n    \ndef getTopLeftCorners(dims, tile_size:int):\n    ''' Make a map of the tiles' locations '''\n    \n    cols, rows = np.ceil(dims/tile_size).astype(int)\n\n    ## M\n    top_left_corners = []\n    for i in range(cols):\n        for j in range(rows):\n            top_left_corners.append( (i*tile_size, j*tile_size) )\n            \n    return np.array(top_left_corners)\n\ndef makeTissueMask(img):\n    ''' Makes a tissue mask. Also filters the green / blue pen markings '''\n\n    hsv = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)\n    \n    # Filter green/blue ink marker annotations\n    lower_red = np.array([120,20,180])\n    upper_red = np.array([190,220,255])\n    tissue_mask = cv2.inRange(hsv, lower_red, upper_red)\n    \n    # Post-process\n    tissue_mask = cv2.dilate(tissue_mask, None, iterations=2)\n    tissue_mask = cv2.morphologyEx(tissue_mask, cv2.MORPH_CLOSE, None)\n    tissue_mask = cv2.medianBlur(tissue_mask, 21)\n    \n    return tissue_mask\n\ndef fillTissueMask(tissue_mask):\n    ''' Fills the holes in a tissue mask by filling the corresponding contours '''\n\n    new_tm = tissue_mask.copy()\n    contours,_ = cv2.findContours(new_tm,0,cv2.CHAIN_APPROX_SIMPLE)\n    cv2.drawContours(new_tm,contours,-1,255,-1)\n    \n    return new_tm\n\ndef estimateWithCircles(arch, radius:int=64):\n\n    arch = arch.squeeze()\n    circle_centers=[arch[0]]\n    for point in arch:\n        if np.all(np.linalg.norm(point-np.array(circle_centers),axis=-1)>radius):\n            circle_centers.append(point)\n    return np.array(circle_centers)\n\ndef getTissuePercentages(tissue_mask, level_offset:int, tile_size:int, top_left_corners):\n    ''' Calculate the tissue percentage per each tile. Tissue ~ non 255 on the red-channel '''\n\n    ds_rate = 4**level_offset\n    tissue_pcts = np.array([tissue_mask[j//ds_rate:(j+tile_size)//ds_rate,\n                                        i//ds_rate:(i+tile_size)//ds_rate].sum()\n                            for (i,j) in top_left_corners])/(255*(tile_size/ds_rate)**2)\n    \n    return tissue_pcts\n\ndef padIfNeeded(img, tgt_width:int=128, tgt_height:int=128, border_color:int=255):\n    ''' Pad images that need padding (padding on right and bottom)'''\n    \n    h,w = img.shape[0:2]\n    \n    if w<tgt_width or h<tgt_height:\n        padded = np.ones((tgt_height,tgt_width,3), dtype='uint8')*border_color\n        padded[:h,:w] = img\n        return padded\n\n    return img\n\ndef distributeIntToChunks(available:int, weights):\n    '''\n        To distribute `available` \"seats\" based on weights, see\n        https://stackoverflow.com/questions/9088403/distributing-integers-using-weights-how-to-calculate\n    '''\n    \n    distributed_amounts = []\n    total_weights = sum(weights)\n    for weight in weights:\n        weight = float(weight)\n        p = weight / total_weights\n        distributed_amount = round(p * available)\n        distributed_amounts.append(distributed_amount)\n        total_weights -= weight\n        available -= distributed_amount\n    return np.int32(distributed_amounts)\n\n\n# ~~~~~~~~~~~~ Slide class ~~~~~~~~~~~~\nclass Slide:\n    \n    def __init__(self,slide_fn, level=2, tile_size=128, mask_fn=None, data_provider=None):\n        self.slide_fn = slide_fn\n        self.level = level\n        self.tile_size = tile_size\n\n        self.img = MultiImage(self.slide_fn)[self.level].copy()\n        self.dims = np.array(self.img.shape[:2][::-1])\n        self.ds_img = MultiImage(self.slide_fn)[2].copy()\n        self.tissue_mask = makeTissueMask(self.ds_img)\n\n\n        self.mask_fn = mask_fn\n        self.data_provider = data_provider\n        if not (self.mask_fn==None or self.data_provider==None):\n            self.mask = MultiImage(mask_fn)[level].sum(axis=-1)\n\n        self.tile_coords = None\n\n\n    def getTileCoords(self, num_tiles:int, sampling_method='skeleton', tissue_th:tuple=(0.2, 0.7), seed=None, ):\n        ''' \n            Find `num_indices` indices with maximal amount of tissue. `tissue_th`  is the slice of tissue percentage ~(min, max)\n            within which we allow the search: sometimes the `max` might not give enough tiles \n            for the mosaic, so we can decrease it gradually until `min` is reached. \n        '''\n\n        assert sampling_method in {'skeleton', 'tissue_pct', 'slic'}\n        \"`sampling_method` should be one of 'skeleton', 'tissue_pct' or 'slic'\"\n        \n        level_offset = 2 - self.level\n        \n        if sampling_method == 'skeleton':\n\n            # Determine skeleton from filled tissue mask\n            filled_tissue_mask = fillTissueMask(self.tissue_mask.copy())\n            filled_tissue_mask = filled_tissue_mask // 255  # needs to be an array of 0s and 1s\n            skeleton = skeletonize(filled_tissue_mask, method='lee')\n            self.skeleton = np.uint8(np.where(skeleton != 0, 255, 0))\n\n            skeleton = cv2.dilate(self.skeleton, None)\n            contours, _ = cv2.findContours(skeleton, 0, cv2.CHAIN_APPROX_NONE)\n\n            # Filter contours based on length\n            arch_lens = []\n            valid_indices = []\n            radius=int( self.tile_size / (4 ** level_offset) )\n            for idx, arch in enumerate(contours):\n                c = cv2.arcLength(arch, False)\n                c = c / 2\n                if c < radius / 3:\n                    continue\n\n                valid_indices.append(idx)\n                arch_lens.append(c)\n\n            if not np.array(valid_indices).size == 0:\n                contours = np.array(contours)[valid_indices]\n\n            # Extract points from the accepted contours\n            weights = np.array(arch_lens) / np.sum(arch_lens)\n            points_per_arch = distributeIntToChunks(num_tiles, weights)  # <- number of points to be extracted\n            for idx, arch in enumerate(contours):\n\n                num_indices = points_per_arch[idx]\n                output = np.zeros_like(skeleton)\n                cv2.drawContours(output, [arch], -1, 1, 1)\n\n                y_, x_ = np.where(output)\n\n                # Simplify the shape by fitting circles\n                arch = np.dstack([x_, y_])\n                cps = estimateWithCircles(arch, radius)\n                cx, cy = cps[..., 0], cps[..., 1]  #\n\n                # Randomly select indices, in case too many; seed if needed\n                if len(cx) > num_indices:\n\n                    # Seed if needed\n                    if not seed == None:\n                        np.random.seed(seed)\n                    indices = sorted(np.random.choice(len(cx), points_per_arch[idx], replace=False))\n                    np.random.seed(None)  # Return clock seed\n\n                    cx, cy = cx[indices], cy[indices]\n\n                # Append to returnables\n                if idx == 0:\n                    intermediate_coords = np.dstack([cx, cy])\n                else:\n                    intermediate_coords = np.hstack([intermediate_coords, np.dstack([cx, cy])])\n\n            # To top-left-corner format\n            final_coords = intermediate_coords.squeeze()*4**level_offset\n            final_coords = final_coords - np.array( [self.tile_size//2,self.tile_size//2])\n\n            self.tile_coords = final_coords.copy()\n\n        elif sampling_method == 'tissue_pct':\n            self.top_left_corners = getTopLeftCorners(self.dims, self.tile_size)\n            self.tissue_pcts = getTissuePercentages(self.tissue_mask, level_offset=level_offset,\n                                                    tile_size=self.tile_size,\n                                                    top_left_corners=self.top_left_corners)\n\n            # Find indices\n            tth_min, tth = tissue_th\n            while len(np.where(self.tissue_pcts > tth)[0]) < num_tiles:\n                if tth <= tth_min:\n                    break\n\n                tth -= 0.05\n\n            # Indices\n            indices = np.where(self.tissue_pcts > tth)[0]\n\n            # Randomly select indices, in case too many; seed if needed\n            if len(indices)>num_tiles:\n                if not seed == None:\n                    np.random.seed(seed)\n                indices = sorted(np.random.choice(indices, num_tiles, replace=False))\n                np.random.seed(None)  # Return clock seed\n\n            self.tile_coords = self.top_left_corners[indices].copy()\n\n        elif sampling_method == 'slic':\n\n            # Determine SLIC clusters from filled tissue mask\n            filled_tissue_mask = fillTissueMask(self.tissue_mask.copy())\n            filled_tissue_mask = filled_tissue_mask // 255  # needs to be an array of 0s and 1s\n\n            segments = seg.slic(self.ds_img, compactness=10, seed_type='nplace', mask=filled_tissue_mask, n_segments=num_tiles,\n                                multichannel=True, recompute_seeds=True, enforce_connectivity=True)\n            indices = [k for k in np.unique(segments) if not k==-1]\n\n            # Randomly select indices, in case too many; seed if needed\n            if len(indices)>num_tiles:\n                if not seed == None:\n                    np.random.seed(seed)\n                indices = sorted(np.random.choice(indices, num_tiles, replace=False))\n                np.random.seed(None)  # Return clock seed\n\n            for i in indices:\n                contours, _ = cv2.findContours(np.uint8(np.where(segments==i, 255, 0)), 0, 1)\n                contours = sorted(contours, key=lambda x: cv2.contourArea(x))[::-1]\n                \n                M = cv2.moments(contours[0])\n                cx = np.int32(M['m10'] / M['m00'])\n                cy = np.int32(M['m01'] / M['m00'])\n\n                if i==0:\n                    intermediate_coords = np.dstack([cx, cy])\n                else:\n                    intermediate_coords = np.hstack([intermediate_coords, np.dstack([cx, cy])]) \n                           \n            # Append more cluster contours if num_tiles has not been reached\n            if len(indices)<num_tiles:\n                enough_tiles = False\n                for i in indices:\n                    contours, _ = cv2.findContours(np.uint8(np.where(segments==i, 255, 0)), 0, 1)\n                    contours = sorted(contours, key=lambda x: cv2.contourArea(x))[::-1]\n                    \n                    for j, cnt in enumerate(contours):\n                        # accept a slic cluster contour if it's area is at least 5% of tile area\n                        if j!=0 and cv2.contourArea(cnt)>(0.05 * self.tile_size / (4**level_offset))**2:\n                            M = cv2.moments(cnt)\n                            cx = np.int32(M['m10'] / M['m00'])\n                            cy = np.int32(M['m01'] / M['m00'])\n\n                            intermediate_coords = np.hstack([intermediate_coords, np.dstack([cx, cy])])\n                            \n                            if intermediate_coords.shape[1] == 12:\n                                enough_tiles = True\n                                break\n                    if enough_tiles:\n                        break\n                            \n            # To top-left-corner format\n            final_coords = intermediate_coords.squeeze() * 4 ** level_offset\n            final_coords = final_coords - np.array([self.tile_size // 2, self.tile_size // 2])\n            self.tile_coords = final_coords.copy()\n\n    def getTiles(self, stack:bool=False, sampling_method:str='skeleton', mosaic_grid:tuple=(4,3), output_tile_size:int=128, tissue_th:tuple=(0.1, 0.7), seed:int=None):\n        ''' Get tiles from the slide and stack into mosaic if needed '''\n\n        # Solve indices to be used in mosaic\n        m, n = mosaic_grid\n        self.getTileCoords(num_tiles=n*m, sampling_method=sampling_method, tissue_th=tissue_th, seed=seed)\n\n        # Read regions\n        output_img = np.ones([n * m, output_tile_size, output_tile_size, 3], dtype='uint8') * 255\n\n        for idx, coord in enumerate(self.tile_coords):\n            x, y = coord\n            left, right = np.int32(np.clip([x, x + self.tile_size], 0, self.dims[0]))\n            bottom, top = np.int32(np.clip([y, y + self.tile_size], 0, self.dims[1]))\n\n            tile = self.img[bottom:top, left:right].copy()\n            tile = padIfNeeded(tile, tgt_width=self.tile_size, tgt_height=self.tile_size)\n            tile = cv2.resize(tile, (output_tile_size,) * 2)\n\n            output_img[idx] = np.uint8(tile)\n\n        if len(self.tile_coords)<m*n:\n            warnings.warn(\"Could not find enough unique tiles for the slide\"\n                          \"(tiles: %s/%s, slide: %s\" %(len(self.tile_coords), m*n, self.slide_fn))\n\n        # Stack to single array of (m,n) tiles if needed\n        if stack:\n            output_img = [np.hstack([output_img[i * n + j] for j in range(n)]) for i in range(m)]\n            output_img = np.vstack(output_img)\n\n        return np.array(output_img)\n\n    def getTilesCancerStatus(self, stack:bool=False, mosaic_grid:tuple=(4,3)):\n        ''' Get tiles from the slide and stack into mosaic if needed '''\n\n        # Solve indices to be used in mosaic\n        m, n = mosaic_grid\n\n        # Read regions\n        output_img = np.zeros([n * m], dtype='uint8')\n\n        for idx, coord in enumerate(self.tile_coords):\n            x, y = coord\n            left, right = np.int32(np.clip([x, x + self.tile_size], 0, self.dims[0]))\n            bottom, top = np.int32(np.clip([y, y + self.tile_size], 0, self.dims[1]))\n\n            tile = self.mask[bottom:top, left:right].copy()\n            tile_cat = tileClassification(tile, self.data_provider)\n\n            output_img[idx] = tile_cat\n\n        # Stack to single array of (m,n) tiles if needed\n        if stack:\n            output_img = [np.hstack([output_img[i * n + j] for j in range(n)]) for i in range(m)]\n            output_img = np.vstack(output_img)\n\n        return np.array(output_img)\n\n    def visualizeCoverage(self, figsize=(16,16)):\n        ''' Visualize the coverage of indices on a slide '''\n\n        background = self.ds_img.copy()\n        foreground = background.copy()\n\n        level_offset = 2-self.level\n\n        for coord in self.tile_coords:\n            x, y = coord\n            left, right = np.int32(np.clip([x, x + self.tile_size], 0, self.dims[0])/(4**level_offset))\n            bottom, top = np.int32(np.clip([y, y + self.tile_size], 0, self.dims[1])/(4**level_offset))\n\n            foreground[bottom:top, left:right] = (0, 255, 0)\n\n        ## Visualize\n        output = cv2.addWeighted(background, 0.7, foreground, 0.3, 0)\n        \n        plt.figure(figsize=figsize)\n        plt.imshow(output)\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class PandaDataset(Dataset):\n    def __init__(self, path, test, level, \n                 tile_size, sz, mosaic_grid,\n                 mean,std,\n                 tissue_th=(0.2,0.7), sampling_method='skeleton', seed=123):\n        self.path = path\n        self.names = list(pd.read_csv(test).image_id)\n        self.level = level\n        self.sz = sz\n        self.tile_size = tile_size\n        self.mosaic_grid = mosaic_grid\n        self.mean = mean\n        self.std = std\n        self.tissue_th = tissue_th\n        self.sampling_method = sampling_method\n        self.seed = seed\n        self.N = mosaic_grid[0]*mosaic_grid[1]\n    def __len__(self):\n        return len(self.names)\n    def __getitem__(self, idx):\n        name = self.names[idx]\n        fn = os.path.join(self.path,name+'.tiff')\n        try:\n            slide = Slide(fn, level=self.level, tile_size=self.tile_size)\n\n            tiles = slide.getTiles(mosaic_grid=self.mosaic_grid,\n                                   stack=False,\n                                   sampling_method=self.sampling_method,\n                                   output_tile_size=self.sz,\n                                   tissue_th=self.tissue_th, seed=self.seed)\n        except Exception as e:\n            print(\"Error in tile generation: %s\" %e)\n            tiles = np.zeros((self.N,self.sz,self.sz,3)).astype(np.float32)\n        tiles = tiles.reshape(-1,self.sz,self.sz,3)\n        tiles = torch.Tensor(1.0 - tiles/255.0)\n        tiles = (tiles - self.mean)/self.std\n        return tiles.permute(0,3,1,2), name","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Slide sampling + Dataloader unit test","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAIN_DATA = '../input/prostate-cancer-grade-assessment/train_images/'\nTRAIN = '../input/prostate-cancer-grade-assessment/train.csv'\ndls = []\ndef to_one(data, sz, mean, std, N):\n    img = torch.stack(data,1)\n    img = img.view(3,-1,N,sz,sz).permute(0,1,3,2,4).contiguous().view(3,-1,sz*N)\n    img = 1.0 - (mean[...,None,None]+img*std[...,None,None])\n    return Image(img)\n\n# don't crash if training data is missing\nif os.path.exists(TRAIN_DATA):\n    \n    for i, _ in enumerate(models):\n        # load training data and correct labels\n        ds = PandaDataset(\n            TRAIN_DATA,\n            TRAIN,\n            level=level[i],\n            tile_size=tile_size[i],\n            sz=sz[i],\n            mosaic_grid=mosaic_grid[i],\n            mean=mean[i],\n            std=std[i],\n            tissue_th=tissue_th,\n            sampling_method=sampling_method[i],\n            seed=seed\n        )\n        dl = DataLoader(ds, batch_size=bs, num_workers=nworkers, shuffle=False)\n        dls.append(dl)\n        train_df = pd.read_csv(TRAIN).set_index('image_id')\n        for (x,y) in dl:\n            break\n        # display batch tiles\n        for img in x:\n            display(to_one(list(img), sz[i], mean[i], std[i], N[i]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model unit test","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# quadratic weights\nw = np.zeros((6,6))\nfor i in range(len(w)):\n    for j in range(len(w)):\n        w[i][j] = float(((i-j)**2)/16)\nw","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def ordinalRegs2cat(arr, classes=[0,1,2,3,4,5]):\n    #mask = arr == 0\n    #return np.clip(np.where(mask.any(1), mask.argmax(1), len(classes)) - 1, classes[0], classes[-1])\n    return np.sum(arr,1)\n\ndef regr2cat(p):\n    for idx,pred in enumerate(p):\n        if   pred < coef[0]: p[idx] = 0\n        elif pred < coef[1]: p[idx] = 1\n        elif pred < coef[2]: p[idx] = 2\n        elif pred < coef[3]: p[idx] = 3\n        elif pred < coef[4]: p[idx] = 4\n        else:                p[idx] = 5\n    return p","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nfrom sklearn.metrics import cohen_kappa_score\n\nRUN_N_TRAIN_SAMPLES = 50\n\n# don't crash if training data is missing\nif os.path.exists(TRAIN_DATA):\n    names,preds = [],[]\n    \n    with torch.no_grad():\n        for i, samples in tqdm(enumerate(zip(*dls))):\n            meta_model_inputs = np.zeros((bs,len(dls)))\n            len_models = 0\n            for j, (x,y) in enumerate(samples):\n                x = x.cuda().float()\n                #dihedral TTA\n                x = torch.stack([x,x.flip(-1),x.flip(-2),x.flip(-1,-2),\n                  x.transpose(-1,-2),x.transpose(-1,-2).flip(-1),\n                  x.transpose(-1,-2).flip(-2),x.transpose(-1,-2).flip(-1,-2)],1)\n                x = x.view(-1,N[j],3,sz[j],sz[j])\n                if is_ordinal[j]:\n                    p = [tensor(ordinalRegs2cat((torch.sigmoid(model(x)) > 0.5).cpu().numpy())).float().cuda() for model in models[j]]\n                elif regr[j]:\n                    p = [regr2cat(model(x)) for model in models[j]]\n                else:\n                    p = [model(x).argmax(-1).float() for model in models[j]]\n                \n                # TTA averages - keep models separate\n                p = torch.stack(p,1).view(bs,len(models[j]), 8, -1) \n                p = p.mean(axis=2)\n                for idx,preds_item in enumerate(p):\n                    softmax_sum = np.zeros((6), np.float) # number of classes\n                    for preds_fold, cm in zip(preds_item, cms[j]):\n                        # get class integer\n                        if not regr and not is_ordinal:\n                            pred_cl = int(preds_fold.argmax(-1))\n                        else:\n                            if   preds_fold < coef[0]: pred_cl = 0\n                            elif preds_fold < coef[1]: pred_cl = 1\n                            elif preds_fold < coef[2]: pred_cl = 2\n                            elif preds_fold < coef[3]: pred_cl = 3\n                            elif preds_fold < coef[4]: pred_cl = 4\n                            else:                pred_cl = 5\n                \n                        cl_probs = cm[:,pred_cl]\n                        softmax_probs = cl_probs / np.sum(cl_probs)\n                        softmax_sum += softmax_probs\n                    meta_model_inputs[idx,j] = np.argmax(softmax_sum)\n                del x\n            \n            ps = np.round(np.average(meta_model_inputs.reshape(bs,-1), axis=1, weights=model_weights))\n            \n            names.append(y)\n            preds.append(torch.tensor(ps))\n            # stop unit testing after enough samples has been collected\n            if i >= (RUN_N_TRAIN_SAMPLES//bs - 1):\n                print(\"Unit test succeeded.\")\n                break\n    gt_grades = np.array([train_df.loc[pred_ids]['isup_grade'].values for pred_ids in names]).flatten()\n    pred_grades = regr2cat(np.array([pred.numpy() for pred in preds]).flatten())\n    qwk = cohen_kappa_score(gt_grades, pred_grades, weights=\"quadratic\")\n    print(\"QWK: {0:.4f}\".format(qwk))\n    assert qwk>0.7, \"QWK test gave low scores. Something is probably wrong.\"\n    \n    \n    names = np.concatenate(names)\n    preds = regr2cat(np.array([pred.numpy() for pred in preds]).flatten()).astype(np.int32)\n    sub_df = pd.DataFrame({'image_id': names, 'isup_grade': preds})\n    print(sub_df.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prediction","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df = pd.read_csv(SAMPLE)\nif os.path.exists(DATA):\n    dls = []\n    \n    for i, _ in enumerate(models):\n        # load training data and correct labels\n        ds = PandaDataset(\n            DATA,\n            TEST,\n            level=level[i],\n            tile_size=tile_size[i],\n            sz=sz[i],\n            mosaic_grid=mosaic_grid[i],\n            mean=mean[i],\n            std=std[i],\n            tissue_th=tissue_th,\n            sampling_method=sampling_method[i],\n            seed=seed\n        )\n        dl = DataLoader(ds, batch_size=bs, num_workers=nworkers, shuffle=False)\n        dls.append(dl)\n    \n    names,preds = [],[]\n\n    with torch.no_grad():\n        for i, samples in tqdm(enumerate(zip(*dls))):\n            meta_model_inputs = np.zeros((bs,len(dls)))\n            len_models = 0\n            for j, (x,y) in enumerate(samples):\n                x = x.cuda().float()\n                #dihedral TTA\n                x = torch.stack([x,x.flip(-1),x.flip(-2),x.flip(-1,-2),\n                  x.transpose(-1,-2),x.transpose(-1,-2).flip(-1),\n                  x.transpose(-1,-2).flip(-2),x.transpose(-1,-2).flip(-1,-2)],1)\n                x = x.view(-1,N[j],3,sz[j],sz[j])\n                if is_ordinal[j]:\n                    p = [tensor(ordinalRegs2cat((torch.sigmoid(model(x)) > 0.5).cpu().numpy())).float().cuda() for model in models[j]]\n                elif regr[j]:\n                    p = [regr2cat(model(x)) for model in models[j]]\n                else:\n                    p = [model(x).argmax(-1).float() for model in models[j]]\n                \n                # TTA averages - keep models separate\n                p = torch.stack(p,1).view(bs,len(models[j]), 8, -1) \n                p = p.mean(axis=2)\n                for idx,preds_item in enumerate(p):\n                    softmax_sum = np.zeros((6), np.float) # number of classes\n                    for preds_fold, cm in zip(preds_item, cms[j]):\n                        # get class integer\n                        if not regr and not is_ordinal:\n                            pred_cl = int(preds_fold.argmax(-1))\n                        else:\n                            if   preds_fold < coef[0]: pred_cl = 0\n                            elif preds_fold < coef[1]: pred_cl = 1\n                            elif preds_fold < coef[2]: pred_cl = 2\n                            elif preds_fold < coef[3]: pred_cl = 3\n                            elif preds_fold < coef[4]: pred_cl = 4\n                            else:                pred_cl = 5\n                \n                        cl_probs = cm[:,pred_cl]\n                        softmax_probs = cl_probs / np.sum(cl_probs)\n                        softmax_sum += softmax_probs\n                    meta_model_inputs[idx,j] = np.argmax(softmax_sum)\n                del x\n            \n            ps = np.round(np.average(meta_model_inputs.reshape(bs,-1), axis=1, weights=model_weights))\n            \n            names.append(y)\n            preds.append(torch.tensor(ps))\n    \n    names = np.concatenate(names)\n    preds = regr2cat(np.array([pred.numpy() for pred in preds]).flatten()).astype(np.int32)\n    sub_df = pd.DataFrame({'image_id': names, 'isup_grade': preds})\n    sub_df.to_csv('submission.csv', index=False)\n    sub_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"sub_df.to_csv(\"submission.csv\", index=False)\nsub_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}