Merge branch 'mac_num_op_testing' into develop

potential fix for MacOS numerical discrepancy

Author: whitews

.github/workflows/build-test-base.yml

@@ -7,7 +7,8 @@ jobs:
   test:
     strategy:
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
+        # 3.13.5 specified due to Windows python linkage issue
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13.5"]
         platform: [ubuntu-latest, macos-latest, windows-latest]
 
     runs-on: ${{matrix.platform}}

.github/workflows/wheels.yml

@@ -4,11 +4,23 @@ on: [workflow_dispatch]
 
 jobs:
   build_wheels:
-    name: Build wheels on ${{ matrix.os }}
-    runs-on: ${{ matrix.os }}
+    name: Build wheels for ${{ matrix.os }}
+    runs-on: ${{ matrix.runs-on }}
     strategy:
       matrix:
-        os: [ubuntu-22.04, windows-2022, macos-12]
+        include:
+          - os: linux-intel
+            runs-on: ubuntu-latest
+          - os: linux-arm
+            runs-on: ubuntu-24.04-arm
+          - os: windows-intel
+            runs-on: windows-latest
+          - os: macos-intel
+            # macos-13 was the last x86_64 runner
+            runs-on: macos-13
+          - os: macos-arm
+            # macos-14+ (including latest) are ARM64 runners
+            runs-on: macos-latest
 
     steps:
       - name: Checkout
@@ -22,13 +34,13 @@ jobs:
           python-version: 3.11
 
       - name: Install cibuildwheel
-        run: python -m pip install cibuildwheel==2.15.0
+        run: python -m pip install cibuildwheel==2.22.0
 
       - name: Build wheels
         env:
           CIBW_SKIP: pp* *-musllinux_* *-manylinux_i686 *-win32  # skip PyPy, musllinux, 32-bit Linux & win32 builds
-          CIBW_ARCHS_MACOS: x86_64 arm64
-          CIBW_BUILD: cp37-* cp38-* cp39-* cp310-* cp311-* cp312-*
+          CIBW_ARCHS_MACOS: auto
+          CIBW_BUILD: cp39-* cp310-* cp311-* cp312-* cp313-*
         run: python -m cibuildwheel --output-dir wheelhouse
 
       - name: Store artifacts

.gitignore

@@ -1,6 +1,9 @@
 # User-specific stuff:
 .idea
 
+# OS-specific files
+.DS_Store
+
 # Compiled Python & C extensions
 *.pyc
 *.so

README.md

@@ -17,7 +17,7 @@ transformation, and gating tasks commonly used within the flow community.
 FlowUtils is part of a suite of Python libraries for analyzing flow 
 cytometry data.  It was developed as an extension to the light-weight 
 [FlowIO library](https://github.com/whitews/FlowIO). FlowIO reads and 
-writes Flow Cytometry Standard (FCS) files, and has zero dependencies. 
+writes Flow Cytometry Standard (FCS) files, and has minimal dependencies. 
 For higher level interaction with flow cytometry data, including GatingML 
 & FlowJo 10 support, see the related 
 [FlowKit project](https://github.com/whitews/FlowKit).

setup.py

@@ -15,15 +15,24 @@
 with open("README.md", "r") as fh:
     long_description = fh.read()
 
-
+# NOTE ON C EXTENSIONS
+# On MacOS, clang 14 and later introduced changes to
+# Fused Multiply-Add (FMA). Prior to version 14, clang
+# was less likely to use FMA. clang 14 began using FMA by
+# default, resulting in numerical differences in some
+# expressions. Specifically, the is_left() function in
+# gate_helpers.c was affected and resulted in failed tests.
+# To enforce consistency across platforms the compile
+# flag '-ffp-contract=off' was added to disable FMA.
+# This fixed the discrepancies on Mac OS.
 logicle_extension = Extension(
     'flowutils.logicle_c',
     sources=[
         'src/flowutils/logicle_c_ext/_logicle.c',
         'src/flowutils/logicle_c_ext/logicle.c'
     ],
     include_dirs=[np.get_include(), 'src/flowutils/logicle_c_ext'],
-    extra_compile_args=['-std=c99', '-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION']
+    extra_compile_args=['-std=c99', '-ffp-contract=off']
 )
 
 gating_extension = Extension(
@@ -33,7 +42,7 @@
         'src/flowutils/gating_c_ext/gate_helpers.c'
     ],
     include_dirs=[np.get_include(), 'src/flowutils/gating_c_ext'],
-    extra_compile_args=['-std=c99', '-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION']
+    extra_compile_args=['-std=c99',  '-ffp-contract=off']
 )
 
 setup(

src/flowutils/_version.py

@@ -1,4 +1,4 @@
 """
 FlowUtils version
 """
-__version__ = "1.2.0b0"
+__version__ = "1.2.2b0"

src/flowutils/gating_c_ext/_gate_helpers.c

@@ -1,5 +1,4 @@
 #include <Python.h>
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION // to avoid a warning
 #include <numpy/arrayobject.h>
 #include "gate_helpers.h"
 

src/flowutils/gating_c_ext/gate_helpers.c

@@ -1,4 +1,5 @@
 #include <stdlib.h>
+#include <stdio.h>
 #include "gate_helpers.h"
 
 double point_is_left(
@@ -9,22 +10,31 @@ double point_is_left(
         double test_point_x,
         double test_point_y
 ) {
-    double is_left = (point_b_x - point_a_x) * (test_point_y - point_a_y) -
-              (test_point_x - point_a_x) * (point_b_y - point_a_y);
+//    double bx_ax_dist = point_b_x - point_a_x;
+//    double tpy_ay_dist = test_point_y - point_a_y;
+//    double tpx_ax_dist = test_point_x - point_a_x;
+//    double by_ay_dist = point_b_y - point_a_y;
+//    double mult1 = bx_ax_dist * tpy_ay_dist;
+//    double mult2 = tpx_ax_dist * by_ay_dist;
+//
+//    double is_left = mult1 - mult2;
+
+    double is_left = ((point_b_x - point_a_x) * (test_point_y - point_a_y)) -
+              ((test_point_x - point_a_x) * (point_b_y - point_a_y));
     return is_left;
 }
 
 int calc_wind_count(double point_x, double point_y, int vert_count, double *poly_vertices) {
-	int wind_count = 0;
-	double vert_a_x;
-	double vert_a_y;
-	double vert_b_x;
-	double vert_b_y;
-	double is_left;
+    int wind_count = 0;
+    double vert_a_x;
+    double vert_a_y;
+    double vert_b_x;
+    double vert_b_y;
+    double is_left;
 
     // loop through all edges of the polygon
     for (int i=0; i<vert_count; i++) {
-        //edge from poly_vertices[i] to poly_vertices[i+1]
+        // edge from poly_vertices[i] to poly_vertices[i+1]
         vert_a_x = poly_vertices[(i * 2) + 0];
         vert_a_y = poly_vertices[(i * 2) + 1];
 
@@ -40,6 +50,8 @@ int calc_wind_count(double point_x, double point_y, int vert_count, double *poly
             if (point_y < vert_b_y) {
                 // point crosses & edge travels upward
                 is_left = point_is_left(vert_a_x, vert_a_y, vert_b_x, vert_b_y, point_x, point_y);
+//                printf("if is_left: %.14f\n", is_left);
+
                 if (is_left > 0) {
                     // point is left of edge
                     wind_count += 1;  // valid 'up' intersection
@@ -49,6 +61,21 @@ int calc_wind_count(double point_x, double point_y, int vert_count, double *poly
             if (vert_b_y <= point_y) {
                 // point crosses & edge travels downward
                 is_left = point_is_left(vert_a_x, vert_a_y, vert_b_x, vert_b_y, point_x, point_y);
+//                printf("vert_a_x: %.15f\n", vert_a_x);
+//                printf("vert_a_y: %.15f\n", vert_a_y);
+//                printf("vert_b_x: %.15f\n", vert_b_x);
+//                printf("vert_b_y: %.15f\n", vert_b_y);
+//                printf("point_x: %.15f\n", point_x);
+//                printf("point_y: %.15f\n", point_y);
+
+                // (point_b_x - point_a_x) * (test_point_y - point_a_y) -
+                // (test_point_x - point_a_x) * (point_b_y - point_a_y)
+//                printf("vert_b_x - vert_a_x: %.15f\n", vert_b_x - vert_a_x);
+//                printf("point_y - vert_a_y: %.15f\n", point_y - vert_a_y);
+//                printf("point_x - vert_a_x: %.15f\n", point_x - vert_a_x);
+//                printf("vert_b_y - vert_a_y: %.15f\n", vert_b_y - vert_a_y);
+//
+//                printf("else is_left: %.15f\n", is_left);
 
                 if (is_left < 0) {
                     // point is right of edge

src/flowutils/gating_py.py

@@ -0,0 +1,118 @@
+"""
+Temporary implementation of points_in_polygon()
+to troubleshoot discrepancies in latest Mac OS
+versions where some test events are incorrectly
+labelled. This will be removed after a solution
+is found but will remain in the repo in the
+'mac_num_op_testing' branch in case the issue
+occurs again in the future.
+"""
+
+def point_is_left(
+        point_a_x,
+        point_a_y,
+        point_b_x,
+        point_b_y,
+        test_point_x,
+        test_point_y
+):
+    is_left = (point_b_x - point_a_x) * (test_point_y - point_a_y) - \
+              (test_point_x - point_a_x) * (point_b_y - point_a_y)
+    return is_left
+
+
+def calc_wind_count(point_x, point_y, vert_count, poly_vertices):
+    wind_count = 0
+
+    # loop through all edges of the polygon
+    for i, (vert_a_x, vert_a_y) in enumerate(poly_vertices):
+        if i >= vert_count - 1:
+            vert_b_x = poly_vertices[0][0]
+            vert_b_y = poly_vertices[0][1]
+        else:
+            vert_b_x = poly_vertices[i + 1][0]
+            vert_b_y = poly_vertices[i + 1][1]
+
+        if vert_a_y <= point_y:
+            if point_y < vert_b_y:
+                # point crosses & edge travels upward
+                is_left = point_is_left(vert_a_x, vert_a_y, vert_b_x, vert_b_y, point_x, point_y)
+                # print("if is_left: %.14f" % is_left)
+
+                if is_left > 0:
+                    # point is left of edge
+                    wind_count += 1  # valid 'up' intersection
+        else:
+            if vert_b_y <= point_y:
+                # point crosses & edge travels downward
+                is_left = point_is_left(vert_a_x, vert_a_y, vert_b_x, vert_b_y, point_x, point_y)
+                # print("vert_a_x: %.14f" % vert_a_x)
+                # print("vert_a_y: %.14f" % vert_a_y)
+                # print("vert_b_x: %.14f" % vert_b_x)
+                # print("vert_b_y: %.14f" % vert_b_y)
+                # print("point_x: %.14f" % point_x)
+                # print("point_y: %.14f" % point_y)
+                # print("else is_left: %.14f" % is_left)
+
+                if is_left < 0:
+                    # point is right of edge
+                    wind_count -= 1  # valid 'down' intersect
+
+    return wind_count
+
+
+def points_in_polygon(poly_vertices, points):
+    """
+    Determines whether points in an array are inside a polygon. Points on the
+    edge of the polygon are considered inclusive. This function uses the
+    winding number method and is robust to complex polygons with crossing
+    boundaries, including the presence of 'holes' created by boundary crosses.
+
+    This implementation is based on the C implementation by Dan Sunday.
+    Original copyright notice:
+        Copyright 2000 softSurfer, 2012 Dan Sunday
+
+    The website containing the above implementation is no longer available,
+    but was archived by the Wayback Machine. The last archived version is
+    available here:
+
+        https://web.archive.org/web/20210504233957/
+
+    :param poly_vertices: Polygon vertices (array of 2-D points)
+    :param points: Points to test for polygon inclusion
+    :return: Array of winding counts for each point. True is inside polygon.
+    """
+    # First, find the polygon's bounding box & store the min/max values
+    min_x = poly_vertices[0][0]
+    max_x = poly_vertices[0][0]
+    min_y = poly_vertices[0][1]
+    max_y = poly_vertices[0][1]
+
+    for i, (vert_x, vert_y) in enumerate(poly_vertices):
+        if vert_x < min_x:
+            min_x = vert_x
+        elif vert_x > max_x:
+            max_x = vert_x
+
+        if vert_y < min_y:
+            min_y = vert_y
+        elif vert_y > max_y:
+            max_y = vert_y
+
+    wind_counts = []
+
+    for i, (point_x, point_y) in enumerate(points):
+
+        if point_x < min_x or point_x > max_x or point_y < min_y or point_y > max_y:
+            wind_count = 0
+        else:
+            wind_count = calc_wind_count(
+                point_x,
+                point_y,
+                len(poly_vertices),
+                poly_vertices
+            )
+
+        wind_counts.append(wind_count)
+
+    return wind_counts

src/flowutils/logicle_c_ext/_logicle.c

@@ -1,5 +1,4 @@
 #include <Python.h>
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION // to avoid a warning
 #include <numpy/arrayobject.h>
 #include "logicle.h"