Merge pull request #379 from jcapriot/numpy2.0_updates

Numpy2.0 updates

Author: jcapriot

.ci/azure/deploy.yml

@@ -1,6 +1,6 @@
 jobs:
 - job:
-  displayName: "Deploy Docs and source"
+  displayName: "Deploy Docs"
   pool:
     vmImage: ubuntu-latest
 
@@ -14,27 +14,13 @@ jobs:
         echo $BRANCH_NAME
       displayName: Report branch parameters
 
-    # Just download all of the items already built
-    - task: DownloadPipelineArtifact@2
-      inputs:
-        buildType: 'current'
-        artifactName: 'wheels'
-        targetPath: 'dist'
-
-    - task: DownloadPipelineArtifact@2
-      inputs:
-        buildType: 'current'
-        artifactName: 'source_dist'
-        targetPath: 'dist'
-
     - task: DownloadPipelineArtifact@2
       inputs:
         buildType: 'current'
         artifactName: 'html_docs'
         targetPath: 'html'
 
     - bash: |
-        ls -l dist
         ls -l html
       displayName: Report downloaded cache contents.
 
@@ -47,14 +33,6 @@ jobs:
         GH_NAME: $(gh.name)
         GH_EMAIL: $(gh.email)
 
-    - bash: |
-        twine upload --skip-existing dist/*
-      displayName: Deploy source and wheels
-      condition: eq(variables.IS_TAG, true)
-      env:
-        TWINE_USERNAME: $(twine.username)
-        TWINE_PASSWORD: $(twine.password)
-
     # upload documentation to discretize-docs gh-pages on tags
     - bash: |
         git clone -q --branch gh-pages --depth 1 https://${GH_TOKEN}@github.com/simpeg/discretize-docs.git

.ci/azure/wheels.yml

@@ -36,27 +36,11 @@ stages:
     jobs:
       - template: .ci/azure/docs.yml
 
-  - stage: BuildWheels
-    dependsOn:
-      - Testing
-      - DocBuild
-    displayName: "Build Wheels"
-    jobs:
-      - template: .ci/azure/wheels.yml
-
-  - stage: BuildSource
+  - stage: Deploy
+    displayName: "Deploy Docs"
     dependsOn:
       - Testing
       - DocBuild
-    displayName: "Build Source distribution"
-    jobs:
-      - template: .ci/azure/sdist.yml
-
-  - stage: Deploy
-    displayName: "Deploy Source, Wheels, and Docs"
-    dependsOn:
-      - BuildWheels
-      - BuildSource
     condition: and(succeeded(), or(eq(variables.IS_TAG, true), eq(variables.IS_MAIN, true)))
     jobs:
       - template: .ci/azure/deploy.yml

azure-pipelines.yml

@@ -1,38 +1,34 @@
 # NumPy include directory
-# The try-except is needed because when things are
-# split across drives on Windows, there is no relative path and an exception
-# gets raised. There may be other such cases, so add a catch-all and switch to
-# an absolute path. Relative paths are needed when for example a virtualenv is
-# placed inside the source tree; Meson rejects absolute paths to places inside
-# the source tree.
-# For cross-compilation it is often not possible to run the Python interpreter
-# in order to retrieve numpy's include directory. It can be specified in the
-# cross file instead:
-#   [properties]
-#   numpy-include-dir = /abspath/to/host-pythons/site-packages/numpy/core/include
-#
-# This uses the path as is, and avoids running the interpreter.
-incdir_numpy = meson.get_external_property('numpy-include-dir', 'not-given')
-if incdir_numpy == 'not-given'
-  incdir_numpy = run_command(py,
-    [
-      '-c',
-      '''import os
+numpy_nodepr_api = ['-DNPY_NO_DEPRECATED_API=NPY_1_22_API_VERSION']
+np_dep = dependency('numpy', required: false)
+if not np_dep.found()
+  # For cross-compilation it is often not possible to run the Python interpreter
+  # in order to retrieve numpy's include directory. It can be specified in the
+  # cross file instead:
+  #   [properties]
+  #   numpy-include-dir = /abspath/to/host-pythons/site-packages/numpy/core/include
+  #
+  # This uses the path as is, and avoids running the interpreter.
+  incdir_numpy = meson.get_external_property('numpy-include-dir', 'not-given')
+  if incdir_numpy == 'not-given'
+    incdir_numpy = run_command(py,
+      [
+        '-c',
+        '''import os
 import numpy as np
 try:
   incdir = os.path.relpath(np.get_include())
 except Exception:
   incdir = np.get_include()
 print(incdir)
-  '''
-    ],
-    check: true
-  ).stdout().strip()
-else
-  _incdir_numpy_abs = incdir_numpy
+    '''
+      ],
+      check: true
+    ).stdout().strip()
+  endif
+  inc_np = include_directories(incdir_numpy)
+  np_dep = declare_dependency(include_directories: inc_np)
 endif
-inc_np = include_directories(incdir_numpy)
-np_dep = declare_dependency(include_directories: inc_np)
 
 # Deal with M_PI & friends; add `use_math_defines` to c_args or cpp_args
 # Cython doesn't always get this right itself (see, e.g., gh-16800), so
@@ -44,7 +40,6 @@ else
   use_math_defines = []
 endif
 
-numpy_nodepr_api = '-DNPY_NO_DEPRECATED_API=NPY_1_9_API_VERSION'
 c_undefined_ok = ['-Wno-maybe-uninitialized']
 
 cython_c_args = [numpy_nodepr_api, use_math_defines]
@@ -54,14 +49,19 @@ if cy_line_trace
   cython_c_args += ['-DCYTHON_TRACE_NOGIL=1']
 endif
 
+cython_args = []
+if cy.version().version_compare('>=3.1.0')
+  cython_args += ['-Xfreethreading_compatible=True']
+endif
+
 cython_cpp_args = cython_c_args
 
 module_path = 'discretize/_extensions'
 
 py.extension_module(
     'interputils_cython',
     'interputils_cython.pyx',
-    include_directories: incdir_numpy,
+    cython_args: cython_args,
     c_args: cython_c_args,
     install: true,
     subdir: module_path,
@@ -71,7 +71,7 @@ py.extension_module(
 py.extension_module(
     'tree_ext',
     ['tree_ext.pyx' , 'tree.cpp', 'geom.cpp'],
-    include_directories: incdir_numpy,
+    cython_args: cython_args,
     cpp_args: cython_cpp_args,
     install: true,
     subdir: module_path,
@@ -82,7 +82,7 @@ py.extension_module(
 py.extension_module(
     'simplex_helpers',
     'simplex_helpers.pyx',
-    include_directories: incdir_numpy,
+    cython_args: cython_args,
     cpp_args: cython_cpp_args,
     install: true,
     subdir: module_path,

discretize/_extensions/meson.build

@@ -13,6 +13,7 @@
     av,
     av_extrap,
     make_boundary_bool,
+    cross2d,
 )
 from discretize.utils.code_utils import deprecate_method, deprecate_property
 
@@ -2237,7 +2238,12 @@ def boundary_edge_vector_integral(self):
         if self.dim > 2:
             Av *= 2
 
-        w_cross_n = np.cross(-w, Av.T @ dA)
+        av_da = Av.T @ dA
+
+        if self.dim == 2:
+            w_cross_n = cross2d(av_da, w)
+        else:
+            w_cross_n = np.cross(av_da, w)
 
         if self.dim == 2:
             return Pe.T @ sp.diags(w_cross_n, format="csr")

discretize/operators/differential_operators.py

@@ -3,7 +3,7 @@
 import numpy as np
 import scipy.sparse as sp
 from scipy.spatial import KDTree
-from discretize.utils import Identity, invert_blocks, spzeros
+from discretize.utils import Identity, invert_blocks, spzeros, cross2d
 from discretize.base import BaseMesh
 from discretize._extensions.simplex_helpers import (
     _build_faces_edges,
@@ -264,7 +264,7 @@ def face_normals(self):  # NOQA D102
             if self.dim == 2:
                 # Take the normal as being the cross product of edge_tangents
                 # and a unit vector in a "3rd" dimension.
-                normal = np.cross(self.edge_tangents, [0, 0, 1])[:, :-1]
+                normal = np.c_[self.edge_tangents[:, 1], -self.edge_tangents[:, 0]]
             else:
                 # define normal as |01 x 02|
                 # therefore clockwise path about the normal is 0->1->2->0
@@ -346,7 +346,7 @@ def edge_curl(self):  # NOQA D102
                 # cp = np.cross(l01, -l20)
                 # cp is a bunch of 1s (where simplices are CCW) and -1s (where simplices are CW)
                 # (but we take the sign here to guard against numerical precision)
-                cp = np.sign(np.cross(l20, l01))
+                cp = np.sign(cross2d(l20, l01))
                 face_areas = face_areas * cp
                 # don't due *= here
 
@@ -812,29 +812,34 @@ def get_interpolation_matrix(  # NOQA D102
                 n_items = self.n_edges
             elif location_type[:-2] == "faces":
                 # grab the barycentric transforms associated with each simplex:
-                ts = transform[
-                    inds,
-                    :,
-                ]
+                ts = transform[inds, :]
                 ts = np.hstack((ts, -ts.sum(axis=1)[:, None]))
                 # use  Whitney 2 - form basis functions for face vector interp
                 faces = self._simplex_faces[inds]
                 areas = self.face_areas
 
                 # [1, 2], [0, 2], [0, 1]
                 if self.dim == 2:
+                    # i j k
+                    # t0 t1 t2
+                    # 0 0 1
+                    # t1 * i - t0 * j
+
                     f0 = (
-                        barys[:, 1] * (np.cross(ts[:, 2], [0, 0, 1])[:, i_dir])
-                        + barys[:, 2] * (np.cross([0, 0, 1], ts[:, 1])[:, i_dir])
-                    ) * areas[faces[:, 0]]
+                        cross2d(barys[:, 1:], ts[:, 1:, 1 - i_dir]) * areas[faces[:, 0]]
+                    )
                     f1 = (
-                        barys[:, 0] * (np.cross(ts[:, 2], [0, 0, 1])[:, i_dir])
-                        + barys[:, 2] * (np.cross([0, 0, 1], ts[:, 0])[:, i_dir])
-                    ) * areas[faces[:, 1]]
+                        cross2d(barys[:, [0, 2]], ts[:, [0, 2], 1 - i_dir])
+                        * areas[faces[:, 1]]
+                    )
                     f2 = (
-                        barys[:, 0] * (np.cross(ts[:, 1], [0, 0, 1])[:, i_dir])
-                        + barys[:, 1] * (np.cross([0, 0, 1], ts[:, 0])[:, i_dir])
-                    ) * areas[faces[:, 2]]
+                        cross2d(barys[:, :-1], ts[:, :-1, 1 - i_dir])
+                        * areas[faces[:, 2]]
+                    )
+                    if i_dir == 1:
+                        f0 *= -1
+                        f1 *= -1
+                        f2 *= -1
                     Aij = np.c_[f0, f1, f2].reshape(-1)
                     ind_ptr = 3 * np.arange(n_loc + 1)
                 else:
@@ -1219,7 +1224,7 @@ def boundary_edge_vector_integral(self):  # NOQA D102
             Pe = self.project_edge_to_boundary_edge
             n_boundary_edges, n_edges = Pe.shape
             index = boundary_face_edges
-            w_cross_n = np.cross(-self.edge_tangents[index], dA)
+            w_cross_n = cross2d(dA, self.edge_tangents[index])
             M_be = (
                 sp.csr_matrix((w_cross_n, (index, index)), shape=(n_edges, n_edges))
                 @ Pe.T

discretize/unstructured_mesh.py

@@ -77,6 +77,7 @@
   invert_blocks
   make_property_tensor
   inverse_property_tensor
+  cross2d
 
 Mesh Utilities
 --------------
@@ -165,6 +166,7 @@
     inv2X2BlockDiagonal,
     makePropertyTensor,
     invPropertyTensor,
+    cross2d,
 )
 from discretize.utils.mesh_utils import (
     meshTensor,