Python: find faces with incorrect normals to flip and flip them
This script finds the average normal from the selected faces, then finds normals that are very different by calculating the dot product between the average and each face normal.

Any face that has a negative dot product will be flipped. This fixes the normals on the faces in your blendfile, although your selection doesn't include all the flipped faces (if selected, they will be fixed).

The script assumes you are in edit mode with all the faces of interest selected.

import bpy, bmesh
from mathutils import Vector

bm = bmesh.from_edit_mesh( bpy.context.object.data )

# Reference selected face indices
bm.faces.ensure_lookup_table()
selFaces = [ f.index for f in bm.faces if f.select ]

# Calculate the average normal vector
avgNormal = Vector()
for i in selFaces: avgNormal += bm.faces[i].normal
avgNormal = avgNormal / len( selFaces )

# Calculate the dot products between the average an each face normal
dots = [ avgNormal.dot( bm.faces[i].normal ) for i in selFaces ]

# Reversed faces have a negative dot product value
reversedFaces = [ i for i, dot in zip( selFaces, dots ) if dot &lt; 0 ]

# Deselect all faces and (later) only select flipped faces as indication of change
for f in bm.faces: f.select = False
bm.select_flush( False )

for i in reversedFaces:
    bm.faces[i].select = True
    bm.faces[i].normal_flip()  # Flip normal

bm.select_flush( True )