我试图通过Cytonize来加速矢量化的numpy操作（或者至少看看我是否可以）。代码计算了给定两个距离矩阵（target_distances和map_distances从某个展平的坐标矢量计算得出）和一些关于距离类型的信息（它的范围可以从0到3，但在本例中我第一次尝试将其全部设为0）的某种压力。numpy版本是www.example.com

import scipy
import numpy as np
from scipy.spatial import distance

def decay(x, s=10):

    return scipy.special.expit(s*x)

def stress(z, target_distances, dim, distance_types, step,
           nrows, ncols):

    row_coords = np.reshape(z[:dim*nrows],(nrows,dim))
    col_coords = np.reshape(z[dim*nrows:dim*(nrows+ncols)],(ncols,dim))

    map_distances = distance.cdist(row_coords, col_coords).copy()
    error = target_distances - map_distances

    I0 = (distance_types==0) | (distance_types==3)
    I2 = distance_types == 2

    return np.sum(error[I0]**2) + np.sum((error[I2] + step)**2*decay(error[I2] + step))

在cythonized版本中，我的pyx文件如下所示（calculus.pyx）

import cython
cimport cython
from libc.stdlib cimport malloc, free

cdef extern from "ctools.h":
    double stress (double*, double**, int**, double, int, int, int)

@cython.boundscheck(False) 
@cython.wraparound(False)
@cython.nonecheck(False)     
def stress_cython(double[:] z, double[:,:] target_distances, int [:,:]
           distance_types, double step, int dim, int nrows, int ncols):

    cdef int i
    cdef int j
    cdef int N = (nrows+ncols)*dim

    z_C = <double*>malloc(sizeof(double)*N)
    target_distances_C = <double**>malloc(sizeof(double*)*nrows)
    distance_types_C = <int**>malloc(sizeof(int*)*nrows)

    for i in range(N):
        z_C[i] = z[i]

    for i in range(nrows):
        target_distances_C[i] = <double *>malloc(sizeof(double)*ncols)
        distance_types_C[i] = <int *>malloc(sizeof(int)*ncols)

        for j in range(ncols):
            target_distances_C[i][j] = target_distances[i,j]
            distance_types_C[i][j] = distance_types[i,j]

    stress_val = stress(z_C, target_distances_C, 
                        distance_types_C, step, nrows, ncols, dim)

    for i in range(nrows):
        free(target_distances_C[i])
        free(distance_types_C[i])

    free(target_distances_C)
    free(distance_types_C)

    return stress_val

外部ctools.c为

# include<stdio.h>

# include<math.h>

double decay(double x){

    return 1/(1+exp(-10*x));

    }

double**dist_pairs(double* z, int nrows, int ncols, int dims)
    {

    int i,j,d;
    double coord1, coord2;
    double**dist;

    dist = (double**)malloc(sizeof(double*)*nrows);
    for (i=0; i<nrows; i++){
        dist[i] = (double *)malloc(sizeof(double)*ncols);
        }

    for (i=0; i<nrows; i++){
       for (j=0; j<ncols; j++){

            dist[i][j] = 0;

            for (d=0; d<dims; d++){

                coord1 = z[d*(nrows+ncols) + i];
                coord2 = z[d*(nrows+ncols) + nrows + j];
                dist[i][j] += pow(coord1-coord2,2.0) ;
                }

            dist[i][j] = sqrt(dist[i][j]);
        }
    }

    return dist;

}

double stress(double* z, double**target_distances, int**distance_types, 
              double step, int nrows, int ncols, int dim){

    int i,j;
    double stress = 0.0;
    double err = 0.0;
    double**map_distances = dist_pairs(z, nrows, ncols, dim);

    for (i=0; i<nrows; i++){
        for (j=0; j<ncols; j++){

            if (distance_types[i][j]==0 || distance_types[i][j]==3){
                stress += pow(target_distances[i][j] - map_distances[i][j],2.0);
                }
            else if (distance_types[i][j]==2){
                err = target_distances[i][j] - map_distances[i][j] + step;
                stress += pow(step,2)*decay(step);

            }

        }

    }

    for (i=0; i<nrows; i++){
        free(map_distances[i]);
    }
    free(map_distances);

    return stress;

}

我做了下面的测试（test.py）

import numpy as np
import time 
from calculus import stress_cython
from pycalculus import stress

nrows = 3000
ncols = 2000
dim = 2
N = 100
is_discrete = 1.0

dt0 = 0
dt1 = 0
difs = []
for i in range(N):
    coordinates = np.random.rand(dim, nrows+ncols)
    coordinates_flat = coordinates.flatten()

    target_distances = np.random.rand(nrows, ncols)
    distance_types = np.zeros((nrows,ncols), dtype='i')

    t0 = time.time()
    stress1 = stress_cython(coordinates_flat, target_distances, distance_types, is_discrete, dim, nrows,
                   ncols)
    t1 = time.time()

    stress2 = stress(coordinates.T.flatten(), target_distances, dim, distance_types, is_discrete,
                     nrows,ncols)

    t2 = time.time()

    dt0 += t1-t0
    dt1 += t2-t1

    difs.append(stress1-stress2)

print(f'cython:{dt0:.2f} python:{dt1:.2f}')

而在我的机器中，计时结果非常相似（~ 4. 47 vs ~ 4. 62）。我不太明白为什么，是不是有一些并行计算在后台进行numpy和scipy？
当我检查我的核心使用情况时，似乎不是这样。我尝试通过

export MKL_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1
export OMP_NUM_THREADS=1

在运行www.example.com之前，在shell中test.py但没有产生任何影响。我在这里缺少的是矢量化操作期间的并行性，还是我只是在特灵优化已经优化的例程，以便在这些库中进行最佳扩展？
ps1：我用参数“-ffast-math”编译了cython代码
ps2：最后，我会让所有内核都处于不同的初始启动条件下，所以如果numpy自动并行化矢量化操作，它对我没有帮助。这就是我试图了解发生了什么的原因之一。