Python update (ARM-software#191)

* Update and tests to release new wrapper on PyPI

Some corrections to cholesky
Tested with Numpy 2.0 (can be use with Numpy >= 1.23.5)

Author: christophe0606

PythonWrapper/build_darwin/create.sh

@@ -3,6 +3,7 @@ cmake -DHOST=YES \
   -DWRAPPER=YES \
   -DCMAKE_POSITION_INDEPENDENT_CODE=YES \
   -DCMAKE_OSX_ARCHITECTURES="arm64;x86_64" \
+  -DCMAKE_OSX_DEPLOYMENT_TARGET="10.9" \
   -DCMSISDSP="path to CMSIS-DSP folder" \
   -DCMAKE_C_FLAGS_RELEASE="-std=c11 -Ofast -ffast-math -DNDEBUG -Wall -Wextra" \
   -DCMAKE_CXX_FLAGS_RELEASE="-fno-rtti -std=c++11 -Ofast -ffast-math -DNDEBUG -Wall -Wextra -Wno-unused-parameter" \

PythonWrapper/cmsisdsp_pkg/src/cmsisdsp_matrix.c

@@ -1908,6 +1908,15 @@ cmsis_arm_mat_cholesky_f32(PyObject *obj, PyObject *args)
     uint32_t row = src_converted.numRows ;
     createf32Matrix(&dst_converted,row,column);
 
+    float32_t *p=dst_converted.pData;
+    for(int r=0;r<row;r++)
+    {
+        for(int c =0; c < column;c++)
+        {
+          *p++=0.0f;
+        }
+    }
+
     arm_status returnValue = arm_mat_cholesky_f32(&src_converted,&dst_converted);
     PyObject* theReturnOBJ=Py_BuildValue("i",returnValue);
     PyObject* dstOBJ=NumpyArrayFromf32Matrix(&dst_converted);
@@ -1939,6 +1948,16 @@ cmsis_arm_mat_cholesky_f64(PyObject *obj, PyObject *args)
     uint32_t row = src_converted.numRows ;
     createf64Matrix(&dst_converted,row,column);
 
+    float64_t *p=dst_converted.pData;
+    for(int r=0;r<row;r++)
+    {
+        for(int c =0; c < column;c++)
+        {
+          *p++=0.0f;
+        }
+    }
+
+
     arm_status returnValue = arm_mat_cholesky_f64(&src_converted,&dst_converted);
     PyObject* theReturnOBJ=Py_BuildValue("i",returnValue);
     PyObject* dstOBJ=NumpyArrayFromf64Matrix(&dst_converted);

PythonWrapper/cmsisdsp_pkg/src/cmsisdsp_module.h

@@ -27,14 +27,15 @@
  */
 #ifndef CMSISMODULE_H
 #define CMSISMODULE_H
-#define NPY_NO_DEPRECATED_API NPY_1_22_API_VERSION
+#define NPY_NO_DEPRECATED_API NPY_1_23_API_VERSION
 
 #include <numpy/numpyconfig.h>
 
-// API version used on google colab
+// Check it is built with right version
+// (should be backward compatible down to 1.23.5)
 // https://github.com/numpy/numpy/blob/main/numpy/_core/include/numpy/numpyconfig.h
-#if (NPY_API_VERSION != 0x0000000F )
-//#error("Error building with wrong NumPy API version")
+#if (NPY_API_VERSION != NPY_2_0_API_VERSION  )
+#error("Error building with wrong NumPy API version")
 #endif
 
 #ifdef WIN

PythonWrapper/examples/debug.py

@@ -1,15 +1,14 @@
-import cmsisdsp as dsp
-import cmsisdsp.fixedpoint as f
-
+import cmsisdsp as dsp 
 import numpy as np
+from numpy.testing import assert_allclose
+from scipy.stats import entropy,tstd, tvar
+from scipy.special import logsumexp
+from scipy.linalg import cholesky,ldl,solve_triangular
 from scipy import signal
-import matplotlib.pyplot as plt
-import scipy.fft
+
 
 import colorama
 from colorama import init,Fore, Back, Style
-from numpy.testing import assert_allclose
-import scipy.spatial.distance as d 
 
 init()
 
@@ -20,61 +19,16 @@ def printSubTitle(s):
     print("\n" + Style.BRIGHT + s + Style.RESET_ALL)
 
 
-def packset(a):
-    b = np.packbits(a)
-    newSize = int(np.ceil(b.shape[0] / 4.0)) * 4
-    c = np.copy(b).astype(np.uint32)
-    c.resize(newSize)
-    #print(c)
-    vecSize = round(newSize/4)
-    c=c.reshape(vecSize,4)
-    #print(c)
-    r = np.zeros(vecSize)
-    result = []
-    for i in range(0,vecSize):
-        print(c[i,:])
-        #print("%X %X %X %X" % (c[i,0],c[i,1],c[i,2],c[i,3]))
-        d = (c[i,0] << 24) | (c[i,1] << 16) | (c[i,2] << 8) | c[i,3] 
-        result.append(np.uint32(d))
-    return(result) 
-
-nb = 34
-#va = np.random.choice([0,1],nb)
-# Array of word32 containing all of our bits
-#pva = packset(va)
-
-
-#vb = np.random.choice([0,1],nb)
-# Array of word32 containing all of our bits
-#pvb = packset(vb)
-#
-va=[1, 0, 1, 0, 1, 1, 1, 0 ,0, 1, 1, 0, 1, 0, 0, 0, 0, 1,0,0,0,1,1,0,1,0,1,0,0,1,1,1,1,1]
-vb=[0,1,0,0,1,0,0,0,1,0,0,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,1,0,0,0,1,0]
-
-va = np.array(va)
-vb = np.array(vb)
+print("Cholesky")
 
-pva=packset(va)
-pvb=packset(vb)
-
-#pva = [np.uint32(167), np.uint32(0)]
-#pvb = [np.uint32(152), np.uint32(0)]
-
-#print(va,pva)
-#print(vb,pvb)
-
-ctt=1.0*np.count_nonzero((va==1) & (vb==1))
-ctf=1.0*np.count_nonzero((va==1) & (vb==0))
-cft=1.0*np.count_nonzero((va==0) & (vb==1))
-
-res=(cft+ctf)/(2*ctt+cft+ctf)
+a=np.array([[4,12,-16],[12,37,-43],[-16,-43,98]])
+ref=cholesky(a,lower=True)
+print(ref)
 
+status,res=dsp.arm_mat_cholesky_f32(a)
 print(res)
+assert_allclose(ref,res,1e-6,1e-6)
 
-
-print("\nDice")
-ref=d.dice(va,vb)
-res=dsp.arm_dice_distance(pva,pvb,nb)
-print(ref)
+status,res=dsp.arm_mat_cholesky_f64(a)
 print(res)
-assert_allclose(ref,res,1e-6)
+assert_allclose(ref,res,1e-10,1e-10)

PythonWrapper/examples/requirements.txt

@@ -0,0 +1,4 @@
+numpy>=1.23.5
+scipy~=1.13.1
+scikit_learn~=1.5.0
+colorama~=0.4.6

PythonWrapper/examples/runall.sh

@@ -0,0 +1,11 @@
+#!/usr/bin/env bash
+
+set -e
+
+python testdistance.py 
+python testdsp6.py 
+python testdsp5.py  
+python example_1_10.py
+python example_1_6.py
+python example_1_9.py
+

PythonWrapper_README.md

@@ -236,6 +236,8 @@ MEL filters are represented as 3 arrays to encode a sparse array.
 
 ## Version 1.9.9:
 * Supports Python 3.12
+* Works with Numpy 2.0
+* Corrections on Cholesky
 
 ## Version 1.9.8:
 * Compute graph API has been removed

Source/MatrixFunctions/arm_mat_cholesky_f16.c

@@ -52,6 +52,10 @@
    * @par
    * If the matrix is ill conditioned or only semi-definite, then it is better using the LDL^t decomposition.
    * The decomposition of A is returning a lower triangular matrix U such that A = L L^t
+   *
+   * @par
+   * The destination matrix should be set to 0 before calling the functions because
+   * the function may not overwrite all output elements.
    */
 
 #if defined(ARM_MATH_MVE_FLOAT16) && !defined(ARM_MATH_AUTOVECTORIZE)

Source/MatrixFunctions/arm_mat_cholesky_f32.c

@@ -60,6 +60,10 @@
    * @par
    * If the matrix is ill conditioned or only semi-definite, then it is better using the LDL^t decomposition.
    * The decomposition of A is returning a lower triangular matrix L such that A = L L^t
+   *
+   * @par
+   * The destination matrix should be set to 0 before calling the functions because
+   * the function may not overwrite all output elements.
    */
 
 #if defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE)

Source/MatrixFunctions/arm_mat_cholesky_f64.c

@@ -50,6 +50,10 @@
    * @par
    * If the matrix is ill conditioned or only semi-definite, then it is better using the LDL^t decomposition.
    * The decomposition of A is returning a lower triangular matrix L such that A = L L^t
+   *
+   * @par
+   * The destination matrix should be set to 0 before calling the functions because
+   * the function may not overwrite all output elements.
    */
 
 #if defined(ARM_MATH_NEON) && !defined(ARM_MATH_AUTOVECTORIZE) && defined(__aarch64__)