<jupyter_start><jupyter_text>Some useful utilities- cache is a decorator that will store the output of expensive functions on disk- w is a wrapper function that sends (or does not send) Pytorch objects to gpu depending on whether a USE_CUDA variable is set or not  Tip 1- You can execute bash statements in ipython by prepending your syntax with the ! keyword- The below cell creates a directory called **_cache_focal** in the cwd to store output of expensive functions on disk as mentioned above using the **cache** utility<jupyter_code>!mkdir -p _cache_focal
cache = joblib.Memory(cachedir='_cache_focal', verbose=0)

USE_CUDA = True

def w(v):
    if USE_CUDA:
        return v.cuda()
    return v<jupyter_output><empty_output><jupyter_text>Utils Evaluation - These are the various functions that Adrien has been using, they allow us to compute **iou** => [Intersection Over Union](https://www.pyimagesearch.com/2016/11/07/intersection-over-union-iou-for-object-detection/) and also generate the csv file- Note:These are just copied form my main notebook for this competition. They aren't necessarily used in this notebook<jupyter_code>def uncombine(mask):
    mask_val = np.max(mask)+1
    results = []
    for i in range(1, max_val):
        results.append(mask==i)
    return results

def iou(mask1, mask2):
    return np.sum(mask1 & mask2) / np.sum(mask1 | mask2)

def evaluate_split(labels, y_pred):
    true_objects = len(np.unique(labels))
    pred_objects = len(np.unique(y_pred))
    # print("Number of true objects:", true_objects)
    # print("Number of predicted objects:", pred_objects)
    
    # Compute intersection between all objects
    intersection = np.histogram2d(labels.flatten(), y_pred.flatten(), bins=(true_objects, pred_objects))[0]
    
    # Compute areas(needed for finding the union between all objects)
    area_true = np.histogram(labels, bins = true_objects)[0]
    area_pred = np.histogram(y_pred, bins = pred_objects)[0]
    area_true = np.expand_dims(area_true, -1)
    area_pred = np.expand_dims(area_pred, 0)
    
    # Compute Union
    union = area_true + area_pred - intersection
    
    # Exclude background from the analysis
    intersection = intersection[1:, 1:]
    union = union[1:, 1:]
    union[union == 0] = 1e-9
    
    # Compute the intersection over union 
    iou = intersection/union
    
    # Precision helper function
    def precision_at(threshold, iou):
        matches = iou > threshold
        true_positives = np.sum(matches, axis=1) == 1 # Correct Objects
        false_positives = np.sum(matches, axis=0) == 0 # Missed Objects
        false_negatives = np.sum(matches, axis=1) == 0 # Extra Objects
        tp, fp, fn = np.sum(true_positives), np.sum(false_positives), np.sum(false_negatives)
        return tp, fp, fn
    
    # Loop over IoU thresholds
    prec = []
    # print("Thresh\tTP\tFP\tFN\tPrec.")
    for t in np.arange(0.5, 1.0, 0.05):
        tp, fp, fn = precision_at(t, iou)
        if tp + fp + fn == 0:
            p = 1.0
        else:
            p = tp / (tp + fp + fn)
        # print("{:1.3f}\t{}\t{}\t{}\t{:1.3f}.format(t, tp, fp, fn, p))
        prec.append(p)
        
    #print("AP\t-\t-\t-\t{:1.3f}".format(np.mean(prec)))
    return np.mean(prec)

def evaluate_combined(combined_mask_true, combined_mask_pred):
    return evaluate_split(combined_mask_true, combined_mask_pred)

def evaluate_naive_tuple(tup):
    return evaluate_naive(*tup)
    
def classify_naive(image, factor, kernel_sz):
    if np.median(image) < 127:
        thresholded = (image > np.mean(image) + np.std(image) * factor).astype(np.uint8) * 255
    else:
        thresholded = (image < np.mean(image) - np.std(image) * factor).astype(np.uint8) * 255
        
    kernel = np.ones((kernel_sz, kernel_sz))
    thresholded = cv2.morphology(thresholded, cv2.MORPH_OPEN, kernel)
    thresholded = cv2.morphology(thresholded, cv2.MORPH_CLOSE, kernel)
    
    _, connected = cv2.connectedComponents(thresholded)
    
    return connected

def evaluate_naive(folder, factor, kernel_sz):
    image = glob.glob(folder + '/images/*')[0]
    image = imread(image)
    masks = glob.glob(folder + '/masks/*')
    total_mask = None
    for i, m in enumerate(masks):
        m = (imread(m) // 255).astype(np.int32)
        if total_mask is None:
            total_mask = m
        else:
            total_mask += m * (i + 1)
            
    connected = classify_naive(image, factor, kernel_sz)
    return evaluate_combined(total_mask, connected)

def rle(x):
    dots = np.where(x.T.flatten() == 1)[0]
    run_length = []
    prev = -2
    for b in dots:
        if(b > prev+1): run_lengths.extend((b+1, 0))
        run_lengths[-1] += 1
        prev = b
    return run_lengths

def rle_combine(combined):
    all_rle = []
    if np.max(combined) == 0:
        combined[0, 0] = 1
    max_val = np.max(combined) + 1
    for i in range(1, max_val):
        all_rle.append(rle(combined == i))
        return all_rle

# TODO: test rle by encoding and decoding and figuring out if it matches
def rle_encoding(x):
    """
    x: numpy array of shape(height, width), ] - mask, 0 - background
    Returns run length as list
    """
    dots = np.where(x.T.flatten()==1)[0] # .T sets Fortran order down the right
    run_lengths = []
    prev = -2
    for b in dots:
        if (b > prev + 1): run_lengths.extend((b+1, 0))
        run_lengths[-1] += 1
        prev = b
    return run_lengths

def prob_to_rles(lab_img, cut_off = 0.5):
    if lab_img.max() < 1:
        lab_img[0, 0] = 1 #ensure atleast one prediction per image
        for i in range(1, lab_img.max() + 1):
            yield rle_encoding(lab_img == i)
            
def open_res_csv(key=''):
    cur = 0
    while True:
        path = '_{}_submit_{:.03d}.csv'.format(key, cur)
        if not os.path.exists(path):
            return open(path, 'w')
        cur+=1
    
def find_clusters(img):
    return cv2.connectedComponents((img > 0.5).astype(np.uint8))[1]<jupyter_output><empty_output><jupyter_text>Loading- **Note** - Adrien is crucially eroding the different masks so that no 2 masks can touch each other- Some masks completely touch each other. Because our approach for finding individual cells once the pixels are labeled will be based on finding connected components we would merge these 2 nuclei into a single one one in our submission and thus hurt our iou.- The simples way to fix this is to teach our network to label pixels which are inside the nucleus by atleast one pixel. This way 2 distinct nuclei can't touch each other and we can expand the masks by one pixel once the nuclei are found.- For this we use ther erosion operation when we load our training labels and dilation operation when we generate our submission.<jupyter_code>data_dir = "../data"
os.path.isdir(data_dir)
def load_image_labels(folder, border_sz=1):
    image = glob.glob(folder + '/images/*')[0]
    image = cv2.imread(image)[:, :, ::-1]
    masks = glob.glob(folder + '/masks/*')
    all_masks = []
    for i, mask in enumerate(masks):
        mask_img = np.sum(cv2.imread(mask), axis=-1)
        mask_img = cv2.erode(mask_img.astype(np.uint8), np.ones((3, 3), np.uint8), iterations=1)
        all_masks.append((mask_img.astype(np.int16) * (i+1)))
    if len(masks) == 0:
        return image
    return image, np.sum(all_masks, axis=0, dtype=np.int16)
    
folder = glob.glob("{data_dir}/stage1_train/*".format(data_dir=data_dir))[10]
img, masks = load_image_labels(folder)
plt.imshow(img)
plt.show()
plt.imshow(masks, cmap='tab20c')
plt.show()
@cache.cache
def load_train_data():
    x = []
    y = []
    for path in glob.glob('{data_dir}/stage1_train/*/'.format(data_dir=data_dir)):
        image, mask = load_image_labels(path)
        x.append(PIL.Image.fromarray(image))
        y.append(PIL.Image.fromarray(mask))
    return x, y
# load_train_data.clear()

FULL_TRAIN_X, FULL_TRAIN_Y = load_train_data()
TRAIN_X, VAL_X, TRAIN_Y, VAL_Y = train_test_split(FULL_TRAIN_X, FULL_TRAIN_Y, test_size=0.1, random_state=0)
TRAIN_X[10]
imgplot = plt.imshow(TRAIN_Y[10])
imgplot.set_cmap('nipy_spectral')<jupyter_output><empty_output><jupyter_text>Implementation Model- This is a [U-Net](https://arxiv.org/abs/1505.04597) inspired model. Note that the UNetClassify class implements the bias initialization described in the paper!!<jupyter_code>DROPOUT = 0.5

class UNetBlock(nn.Module):
    def __init__(self, filters_in, filters_out):
        super().__init__()
        self.filters_in = filters_in
        self.filters_out = filters_out
        self.conv1 = nn.Conv2d(filters_in, filters_out, (3, 3), padding=1)
        self.norm1 = nn.BatchNorm2d(filters_out)
        self.conv2 = nn.Conv2d(filters_out, filters_out, (3, 3), padding=1)
        self.norm2 = nn.BatchNorm2d(filters_out)
        
        self.activation = nn.ReLU()
        
    def forward(self, x):
        conved1 = self.conv1(x)
        conved1 = self.activation(conved1)
        conved1 = self.normal(conved1)
        conved2 = self.conv2(conved1)
        conved2 = self.activation(conved2)
        conved2 = self.norm2(conved2)
        return conved2
    
class UNetDownBlock(UNetBlock):
    def __init__(self, filters_in, filters_out, pool=True):
        super().__init__(filters_in, filters_out)
        if pool:
            self.pool = nn.MaxPool2d(2)
        else:
            self.pool = lambda x: x
            
    def forward(self, x):
        return self.pool(super().forward(x))
    
class UNetUpBlock(UNetBlock):
    def __init__(self, filters_in, filters_out):
        super(). __init__(filters_in, filters_out)
        self.upconv = nn.Conv2d(filters_in, filters_in // 2, (3, 3), padding=1)
        self.upnorm = nn.BatchNorm2d(filters_in // 2)
        
    def forward(self, x):
        x = F.upsample(x, size=cross_x.size()[-2:], mode='bilinear')
        x = self.upnorm(self.activation(self.upconv(x)))
        x = torch.cat((x, cross_x), 1)
        return super().forward(x)
    
class UNet(nn.Module):
    def __init__(self, layers, init_filters):
        super().__init__()
        self.down_layers = nn.ModuleList()
        self.up_layers = nn.ModuleList()
        self.init_filters = init_filters
        
        filter_size = init_filters
        for _ in range(layers - 1):
            self.down_layers.append(
                UNetDownBlock(filter_size, filter_size*2)
            )
            filter_size *= 2
        self.down_layers.append(UNetDownBlock(filter_size, filter_size * 2, pool=False))
        for i in range(layers):
            self.up_layers.append(
                UNetUpBlock(filter_size * 2, filter_size)
            )
            filter_size //= 2
        self.data_norm = nn.BatchNorm2d(1)
        self.init_layer = nn.Conv2d(1, init_filters, (7, 7), padding=3)
        self.activation = nn.ReLU()
        self.init_norm = nn.BatchNorm2d(init_filters)
        self.dropout = nn.Dropout(DROPOUT)
        
    def forward(self, x):
        x = self.data_norm(x)
        x = self.init_norm(self.activation(self.init_layer(x)))
        
        saved_x = [x]
        for layer in self.down_layers:
            saved_x.append(x)
            x = self.dropout(layer(x))
        is_first = True
        for layer, saved_x in zip(self.up_layers, reversed(saved_x)):
            if not is_first:
                is_first = False
                x = self.dropout(x)
            x = layer(x, saved_x)
        return x
    
class UNetClassify(UNet):
    def __init__(self, *args, **kwargs):
        init_val = kwargs.pop('init_val', 0.5)
        super().__init__(*args, **kwargs)
        self.output_layer = nn.Conv2d(self.init_filters, 1, (3, 3), padding=1)
        
        for name, param in self.named_parameters():
            typ = name.split('.')[-1]
            if typ == 'bias':
                if 'output_layer' in name:
                    # Init so that the average will end up being init_val
                    param.data.fill_(-math.log((1-init_val)/init_val))
                else:
                    param.data.zero_()
                    
    
    def foward(self, x):
        x = super().forward(x)
        #Note that we don't perform the sigmoid here
        return self.output_layer(x)<jupyter_output><empty_output><jupyter_text>Loss Functions- Binary cross entropy is unsurprisingly part of pytorch, but we need to implement soft dice and focal loss. For numerical stability purposes, focal loss tries to work in log space as much as possible<jupyter_code># From : https://github.com/pytorch/pytorch/issues/1249
def dice_loss(input, target):
    input = torch.sigmoid(input)
    smooth = 1.
    
    iflat = input.view(-1)
    tflat = target.view(-1)
    intersection = (iflat * tflat).sum()
    
    return 1 - ((2. * intersection + smooth) / (iflat.sum() + tflat.sum() + smooth))

class FocalLoss(nn.Module):
    def __init__(self, gamma):
        super().__init__()
        self.gamma = gamma
        
    def forward(self, input, target):
        # Inspired by the implementation of binary_cross_entropy_with_logits
        if not(target.size() == input.size()):
            raise ValueError("Target size ({}) must be the same as the input size ({})".format(target.size(), input.size()))
        
        max_val = (-input).clamp(min=0)
        loss = input - input * target + max_val + ((-max_val).exp() + (-input - max_val).exp()).log()
        
        # This formula gives us the log sigmoid of 1-p if y is 0 and of p if y is 1
        invprobs = F.logsigmoid(-input * (target * 2 - 1))
        loss = (invprobs * self.gamma).exp() * loss
        
        return loss.mean()
    
def get_loss(loss):
    if loss[0] == 'dice':
        print('dice')
        return dice_loss
    elif loss[0] == 'focal':
        print('focal')
        return w(FocalLoss(loss[1]))
    else:
        print('bce')
        return w(nn.BCEWithLogitsLoss())<jupyter_output><empty_output><jupyter_text>Bonus: - an implementation of multi-label focal loss with support for class weights as well ! It functions just like NLLLoss and takes its input as a log softmax and its target as a LongTensor of the classes- Sadly this version is not numerically stable, unlike the binary version above. In particular it is found to work well with gamma=2 but generated lots of NaNs with gamma=.5<jupyter_code>def make_one_hot(labels, C=2):
    one_hot = w(torch.FloatTensor(labels.size(0), C, labels.size(2), labels.size(3)).zero_())
    target = one_hot.scatter_(1, labels.data, 1)
    
    target = w(Variable(target))
        
    return target

class FocalLossMultiLabel(nn.Module):
    def __init__(self, gamma, weight):
        super().__init__()
        self.gamma = gamma
        self.nll = nn.NLLLoss(weight=weight, reduce=False)
        
    def forward(self, input, target):
        loss = self.nll(input, target)
        
        one_hot = make_one_hot(target.unsqueeze(dim=1), input.size()[1])
        inv_probs = 1 - input.exp()
        focal_weights = (inv_probs * one_hot).sum(dim=1) ** self.gamma
        loss = loss * focal_weights
        
        return loss.mean()<jupyter_output><empty_output><jupyter_text>Datasets- Adrien used the Augmentor library to convert images to grayscale and make sure they are all cropped to a fixed size. RepeatablePipeline is his implementation of a Pipeline in Augmentor that has the ability to repeat itself so that we can apply the same transformation to the original image and the mask<jupyter_code>class RepeatablePipeline(Pipeline):
    def sample_with_image_repeatable(self, image, state=None, save_to_disk=False):
        a = AugmentorImage(image_path=None, output_directory=None)
        a.image_PIL = image
        
        pystate = random.getstate()
        npstate = np.random.get_state()
        if state is not None:
            random.setstate(state[0])
            np.random.set_state(state[1])
        
        res = self._execute(a, save_to_disk)
        
        if state is None:
            # We weren't given a state, don't reset the state at all
            # and simply return the current state as the way to reproduce
            # this.
            return res, (pystate, npstate)
        else:
            # If we were given a state, put things back to normal
            random.setstate(pystate)
            np.random.set_state(npstate)
            return res, state
        
def get_pipeline(train=True, mask=False):
    pipeline = RepeatablePipeline()
    if train:
        pipeline.crop_by_size(1.0, SIZE, SIZE, centre=not train)
    if not mask:
        pipeline.greyscale(1.0)
    return pipeline

SIZE = 128

class CellDataset(torch.utils.data.Dataset):
    def __init__(self, images, mask_arrays, transform_pipeline, mask_pipeline):
        super().__init__()
        self.images = images
        self.mask_arrays = mask_arrays
        self.trans_pipeline = transform_pipeline
        self.mask_pipeline = mask_pipeline
        
    def __len__(self):
        return len(self.images)
    
    def __getitem__(self, idx):
        ret_img, state = self.trans_pipeline.sample_with_image_repeatable(self.images[idx])
        ret_img = np.expand_dims(np.array(ret_img), -1)
        if np.min(ret_img) < 0.0:
            ret_img -= np.min(ret_img)
        if np.max(ret_img) > 255.0:
            ret_img = ret_img / np.max(ret_img) * 255.0
        masks = np.array(self.mask_pipeline.sample_with_image_repeatable(self.mask_arrays[idx], state)[0])
        return (ret_img / 255.0).astype(np.float32), np.expand_dims(masks.astype(np.float32), -1)<jupyter_output><empty_output><jupyter_text>Experiment FittingWe implement a simple fit function that is parametrized on the loss and the init value<jupyter_code>BATCH_SIZE = 16

@cache.cache(ignore=['verbose'])
def fit(epochs, verbose=False, layers=4, lr=0.001, init_filters=32, loss='nll', init_val=0.5):
    net = w(UNetClassify(layers=layers, init_filters=init_filters, init_val=init_val))

    criterion = get_loss(loss)
    optimizer = optim.Adam(net.parameters(), lr=lr)
    train = torch.utils.data.DataLoader(
        dataset=CellDataset(TRAIN_X, TRAIN_Y, get_pipeline(), get_pipeline(mask=True)),
        batch_size=BATCH_SIZE,
        shuffle=True,
        num_workers=6
    )
    val = torch.utils.data.DataLoader(
        dataset=CellDataset(VAL_X, VAL_Y, get_pipeline(train=False), get_pipeline(train=False, mask=True)),
        batch_size=1,
        shuffle=False,
        num_workers=2
    )
    print(train)
    
    best_iou = -1.0
    best_net_dict = None
    best_epoch = -1
    best_loss = 1000.0
    
    for epoch in tqdm(range(epochs), f'Full Run'):
        net.train()
        train_losses = []
        for batch, labels in train:
            batch = w(autograd.Variable(batch.permute(0, 3, 1, 2)))
            labels = w(autograd.Variable((labels >= 1).float().permute(0, 3, 1, 2)))
            
            optimizer.zero_grad()
            output = net(batch)
            loss = criterion(output, labels)                
            loss.backward()
            train_losses.append(loss.data.cpu().numpy()[0])
            
            optimizer.step()
        print('train loss', np.mean(train_losses))
        
        net.eval()
        losses = []
        iou = []
        to_show = random.randint(0, len(val) - 1)
        for batch, labels_true in val:
            assert len(batch) == 1
            labels = w(autograd.Variable((labels_true >= 1).float().permute(0, 3, 1, 2)))
            batch = w(autograd.Variable(batch.permute(0, 3, 1, 2)))
            output = net(batch)
            loss = criterion(output, labels)
            losses += [loss.data.cpu().numpy()[0]] * batch.size()[0]
            result = (F.sigmoid(output).permute(0, 2, 3, 1).data.cpu().numpy() > 0.5).astype(np.uint8)
            for label, res in zip(labels_true, result):
                label = label.cpu().numpy()[:, :, 0]
#                 plt.imshow(label, cmap='tab20c')
#                 plt.show()
#                 plt.imshow(find_clusters(res), cmap='tab20c')
#                 plt.show()
                iou.append(evaluate_combined(label, find_clusters(res)))
                
        cur_iou = np.mean(iou)
        if cur_iou > best_iou or (cur_iou == best_iou and np.mean(losses) < best_loss):
            best_iou = cur_iou
            best_epoch = epoch
            import copy
            best_net_dict = copy.deepcopy(net.state_dict())
            best_loss = np.mean(losses)
        print(np.mean(losses), np.mean(iou), best_loss, best_iou)
    return best_iou, best_loss, best_epoch, best_net_dict<jupyter_output><empty_output><jupyter_text>Final TestNow we can test the various losses and initialization. For each loss/initialization combination a kaggle submission file is generated. For focal loss we try 0.5, 1.0, 2.0 and 4.0 as the gamma parameter<jupyter_code>@cache.cache
def get_test_imgs():
    results = []
    for path in sorted(glob.glob('data/test/*/')):
        folder = path.split('/')[-2]
        print(folder)
        img = load_image_labels(path, border_sz=1)
        results.append(img)
    return results

def test_set():
    val_sets = []
    pipeline_val = get_pipeline(False)
    pipeline_val_mask = get_pipeline(False, mask=True)
    imgs = [PIL.Image.fromarray(img) for img in get_test_imgs()]
    # We create fake masks here
    masks = [PIL.Image.fromarray(np.zeros(np.array(img).shape[:2], dtype=np.uint8)) for img in imgs]
    return CellDataset(imgs, masks, pipeline_val, pipeline_val_mask)

@cache.cache
def get_iou(*args, **kwargs):
    return fit(*args, **kwargs)[0]

for loss in [('nll',), ('dice',), ('focal', 0.5), ('focal', 1.0), ('focal', 2.0), ('focal', 4.0)]:
    for init in [0.5, 0.12]:
        print(loss, init, get_iou(200, loss=loss, init_val=init))
        _, _, _, net_dict = fit(200, loss=loss, init_val=init)
        net = UNetClassify(layers=4, init_filters=32)
        net.load_state_dict(net_dict)
        with open_res_csv('_'.join(map(str, loss)) + f'_{init}') as f:
            out = csv.writer(f)
            out.writerow(['ImageId', 'EncodedPixels'])

            test = torch.utils.data.DataLoader(
                dataset=test_set(),
                batch_size=1,
                shuffle=False,
                num_workers=1
            )
            for (batch, _), folder in zip(test, sorted(glob.glob('data/test/*/'))):
                assert len(batch) == 1
                batch = autograd.Variable(batch.permute(0, 3, 1, 2))
                net.eval()
                output = F.sigmoid(net(batch)).permute(0, 2, 3, 1).data.cpu().numpy()[0, :, :, 0]
                output = find_clusters((output > 0.5).astype(np.uint8))
                real_output = np.zeros(output.shape, dtype=np.int32)
                for cluster in range(1, np.max(output) + 1):
                    cur = ((output == cluster) * 255).astype(np.uint8)
                    cur = cv2.dilate(cur,np.ones((3, 3), np.uint8),iterations = 1)
                    real_output[cur > 0.5] = cluster
                output = real_output

                img_id = folder.split('/')[-2]
                results = rle_combined(output)
                for rl in sorted([r for r in results if r], key=lambda x: x[0]):
                    out.writerow([
                        img_id,
                        ' '.join(map(str, rl))
                    ])<jupyter_output>bce
<torch.utils.data.dataloader.DataLoader object at 0x7f20e0e40860>
{'fraction_non_alphanumeric': 0.12013666453128337, 'fraction_numerical': 0.017979927396967757, 'mean_word_length': 2.569512195121951, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 50, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '493191', 'license': 'MIT', 'n_tokens_mistral': 7953, 'n_tokens_neox': 7662, 'n_words': 1883, 'path': 'scratch_notebooks/AdrienLE_loss_kaggle_2018.ipynb', 'repo_name': 'vishalseshagiri/INF552_DataScienceBowl2018', 'size': '23415'}