Problem: You are given a binary matrix \(M\) of size \((n,m)\) composed only of 0 and 1.
This matrix represents a map where 1 is the land and 0 is the sea. Each 1 (=land) is connected to other 1s (=land) only vertically and horizontally (i.e. not diagonally). In this exercise, \(M\) can be assumed to be a numpy array.
You can switch only one 0 into a 1 so that your map has the biggest island possible given \(M\). Write an efficient algorithm that finds the location of this 0 to switch and returns the size of the biggest island.

The solution described below has a time complexity of \(O(mn)\).

import numpy as np
from itertools import combinations

# Main function
def find_biggest_island(M):
    # First find the islands in the map
    islands = find_islands(M)
    biggest_island = max(len(i) for i in islands.values())
    # Compute the land_to_island dictionary will be 
    # useful later and make our algorithm more efficient
    land_to_island = create_land_to_island(islands)
    # Figure out which 0 will create the biggest island
    # Such a 0 needs to be connected to at least 2 lands
    for i in range(M.shape[0]):
        for j in range(M.shape[1]):
            if M[i,j] == 0:
                # Find how big the new biggest isalnd would be by switching this zero
                new_island_size = potential_new_island(M, i, j, land_to_island, islands)
                if new_island_size > biggest_island:
                    biggest_island = new_island_size
    return biggest_island

def find_islands(M):
    """
    Take the matrix M as input and returns a list of islands
    An island is defined as a list of tuples (i,j) which are the
    list of 1 that are "connected" vertically or horizontally
    Returns a dictionary of int -> List[(i,j)]
    """
    visited = M.copy()
    islands = {}
    island_number = 0
    for i in range(visited.shape[0]):
        for j in range(visited.shape[1]):
            if visited[i,j] == 1:
                island = bfs(M,i,j)
                islands[island_number] = island
                island_number += 1
                # cross out land already visited
                for u,v in island:
                    visited[u,v] = 0
    return islands

def bfs(M, i, j):
    """
    Breadth First Search algorithm implementation on a binary matrix
    start_loc = (i,j)
    Returns an island = list of tuples
    """
    if M[i, j] == 0:
        return []
    visited = set([(i,j)])
    queue = [(i,j)]
    while len(queue) > 0:
        x, y = queue.pop(0)
        neighbors = neighboring_land(M, x, y)
        for neighbor in neighbors:
            if neighbor not in visited:
                visited.add(neighbor)
                queue.append(neighbor)
    return list(visited)

def neighboring_land(M, i, j):
    """
    Returns the list of connected land from a location
    """
    m, n = M.shape
    neighboring_lands = []
    if i > 0 and M[i-1,j] == 1:
        neighboring_lands.append((i-1,j))
    if i < m - 1 and M[i+1,j] == 1:
        neighboring_lands.append((i+1,j))
    if j > 0 and M[i,j-1] == 1:
        neighboring_lands.append((i,j-1))
    if j < n - 1 and M[i,j+1] == 1:
        neighboring_lands.append((i,j+1))
    return neighboring_lands

def create_land_to_island(islands):
    """
    Creates a dictionary of (i,j) -> int (=island size the land belongs to)
    """
    land_to_island = {}
    for i in islands:
        lands = islands[i]
        for land in lands:
            land_to_island[land] = i
    return land_to_island

def potential_new_island(M, i, j, land_to_island, islands):
    neighbors = neighboring_land(M, i, j)
    if len(neighbors) == 0:
        # Not even worth trying
        new_island_size = 0
    elif len(neighbors) == 1:
        neighbor_island = land_to_island[neighbors[0]]
        new_island_size = 1 + len(islands[neighbor_island])
    else:
        new_island_size = 0
        # Take all combintaions of neighbors
        # Have to deal with special case where two neighbors
        # may or may not belong to the same island
        for comb in combinations(neighbors, 2):
            neighbor1 = comb[0]
            neighbor2 = comb[1]
            neighbor_island1 = land_to_island[neighbor1]
            neighbor_island2 = land_to_island[neighbor2]
            if neighbor_island1 != neighbor_island2:
                potential_island_size = 1 + len(islands[neighbor_island1]) + len(islands[neighbor_island2])
            else:
                potential_island_size = 1 + len(islands[neighbor_island1])
            # Take max
            if potential_island_size > new_island_size:
                new_island_size = potential_island_size
    return new_island_size
    
# Write test cases
M = np.array([[0,0,1],
              [1,0,1],
              [1,1,0]])
answer = 6

find_biggest_island(M)