Problem when trying to cut linestring at multiple points after conversion from Multi-linestring

Question

This is my script:

def cut_line_at_points(line, points):
    '''
    line is gdf format
    points are series of points in gdf format
    '''
    # union all geometries - convert the gdf to shpely (multi)line and (multi)point
    line = line.geometry.unary_union
    points = points.geometry.unary_union

    # check if line in multi-part and merge all parts to one linestring
    if line.geom_type != 'LineString':
        # Put the sub-line coordinates into a list of sublists
        outcoords = [list(i.coords) for i in line]
        # Flatten the list of sublists and use it to make a new line
        line = LineString([i for sublist in outcoords for i in sublist])

    # First coords of line
    coords = list(line.coords)

    # Keep list coords where to cut (cut = 1)
    cuts = [0] * len(coords)
    cuts[0] = 1
    cuts[-1] = 1

    # Add the coords from the points
    coords += [list(p.coords)[0] for p in points]
    cuts += [1] * len(points)

    # Calculate the distance along the line for each point
    dists = [line.project(Point(p)) for p in coords]

    # create new df and sort again if geom_type was LineString

    # sort the coords/cuts based on the distances
    coords = [p for (d, p) in sorted(zip(dists, coords))]
    cuts = [p for (d, p) in sorted(zip(dists, cuts))]

    # generate the Lines
    lines = []
    for i in range(len(coords)-1):
        if cuts[i] == 1:
            # find next element in cuts == 1 starting from index i + 1
            try:
                j = cuts.index(1, i + 1)
                lines.append(LineString(coords[i:j+1]))
            except:
                cuts[-1] = 1
                j = cuts.index(1, i + 1)
                lines.append(LineString(coords[i:j+1]))

    # transform Geometry Collection to GeoDataFrame
    segments = [feature for feature in lines]
    # set the geometry and other columns
    gdf_segments = gpd.GeoDataFrame(list(range(len(segments))), geometry=segments)
    gdf_segments.columns = ['index', 'geometry']

    return gdf_segments

everything is working really well when it's a shapely.LineString

but when I'm trying to cut the Multi-linestring after the flattening it all goes south:

I found out that when Calculating the distance along the line it scrambles the points distances, in this line:

dists = [line.project(Point(p)) for p in coords]

but I don't have any idea how to fix it. I use Python 3.7 and Shapely.