def find_boundaries(boundaries, ecs, monomer, multimer, symmetric):
    """
    Identify axis boundaries for contact map plots

    Parameters
    ----------
    boundaries : {"union", "intersection", "ecs", "structure"} or tuple
             or list(tuple, tuple)
        Set axis range (min/max) of contact map as follows:

        * "union": Positions either in ECs or 3D structure
        * "intersection": Positions both in ECs and 3D structure
        * "ecs": Positions in ECs
        * "structure": Positions in 3D structure
        * tuple(float, float): Specify upper/lower bound manually
        * [(float, float), (float, float)]: Specify upper/lower bounds
          for both x-axis (first tuple) and y-axis (second tuple)

    ecs : pandas.DataFrame
        Table of evolutionary couplings to plot (using columns
        "i" and "j")
    monomer : evcouplings.compare.distances.DistanceMap
        Monomer distance map (intra-chain distances)
    multimer : evcouplings.compare.distances.DistanceMap
        Multimer distance map (multimeric inter-chain distances)
    symmetric : bool
        Sets if distance maps and ECs are symmetric (intra-chain or homomultimer),
        or not (inter-chain).

    Returns
    -------
    (min_x, max_x) : (float, float)
        First and last position on x-axis
    (min_y, max_y) : (float, float)
        First and last position on y-axis
    """

    def _find_pos(axis):
        """
        Find first and last index along a single contact map axis
        """
        # determine what sets of positions are for ECs/contact maps
        ec_pos = set()
        monomer_pos = set()
        multimer_pos = set()

        # need to merge i and j here if symmetric
        if ecs is not None:
            if symmetric:
                ec_pos = set(ecs.i.astype(int)).union(ecs.j.astype(int))
            else:
                ec_pos = set(getattr(ecs, axis).astype(int))

        if monomer is not None:
            monomer_pos = set(
                getattr(monomer, "residues_" + axis).id.astype(int)
            )

        if multimer is not None:
            multimer_pos = set(
                getattr(multimer, "residues_" + axis).id.astype(int)
            )

        structure_pos = monomer_pos.union(multimer_pos)

        # maximum ranges spanned by structure or ECs
        # if any of the sets is not given, revert to
        # the other set of positions in else case
        # (in these cases, union and intersection will
        # be trivially the one set that is actually defined)
        if len(ec_pos) > 0:
            min_ec, max_ec = min(ec_pos), max(ec_pos)
        else:
            min_ec, max_ec = min(structure_pos), max(structure_pos)

        if len(structure_pos) > 0:
            min_struct, max_struct = min(structure_pos), max(structure_pos)
        else:
            min_struct, max_struct = min(ec_pos), max(ec_pos)

        # determine and set plot boundaries
        if boundaries == "union":
            min_val = min(min_ec, min_struct)
            max_val = max(max_ec, max_struct)
        elif boundaries == "intersection":
            min_val = max(min_ec, min_struct)
            max_val = min(max_ec, max_struct)
        elif boundaries == "ecs":
            min_val = min_ec
            max_val = max_ec
        elif boundaries == "structure":
            min_val = min_struct
            max_val = max_struct
        else:
            raise ValueError(
                "Not a valid value for boundaries: {}".format(
                    boundaries
                )
            )

        return min_val, max_val

    # check first if range is specified manually
    if isinstance(boundaries, tuple):
        if len(boundaries) != 2:
            raise ValueError(
                "boundaries must be a tuple with 2 elements (min, max)."
            )
        min_x, max_x = boundaries
        min_y, max_y = boundaries
    elif isinstance(boundaries, list):
        if len(boundaries) != 2 or len(boundaries[0]) != 2 or len(boundaries[1]) != 2:
            raise ValueError(
                "boundaries must be a list of 2 tuples with 2 elements "
                "[(min_x, max_x), (min_y, max_y)]."
            )
        min_x, max_x = boundaries[0]
        min_y, max_y = boundaries[1]
    else:
        min_x, max_x = _find_pos("i")
        min_y, max_y = _find_pos("j")

    return (min_x, max_x), (min_y, max_y)