C++ PyArray_DATA函数代码示例

//.........这里部分代码省略.........
        int res = PyArray_DTypeFromObjectHelper(objects[i], maxdims - 1,
                                                out_dtype, string_type);
        if (res < 0) {
            Py_DECREF(seq);
            goto fail;
        }
        else if (res > 0) {
            Py_DECREF(seq);
            return res;
        }
    }

    Py_DECREF(seq);

    return 0;


promote_types:
    /* Set 'out_dtype' if it's NULL */
    if (*out_dtype == NULL) {
        if (!string_type && dtype->type_num == NPY_STRING) {
            Py_DECREF(dtype);
            return RETRY_WITH_STRING;
        }
        if (!string_type && dtype->type_num == NPY_UNICODE) {
            Py_DECREF(dtype);
            return RETRY_WITH_UNICODE;
        }
        *out_dtype = dtype;
        return 0;
    }
    /* Do type promotion with 'out_dtype' */
    else {
        PyArray_Descr *res_dtype = PyArray_PromoteTypes(dtype, *out_dtype);
        Py_DECREF(dtype);
        if (res_dtype == NULL) {
            return -1;
        }
        if (!string_type &&
                res_dtype->type_num == NPY_UNICODE &&
                (*out_dtype)->type_num != NPY_UNICODE) {
            Py_DECREF(res_dtype);
            return RETRY_WITH_UNICODE;
        }
        if (!string_type &&
                res_dtype->type_num == NPY_STRING &&
                (*out_dtype)->type_num != NPY_STRING) {
            Py_DECREF(res_dtype);
            return RETRY_WITH_STRING;
        }
        Py_DECREF(*out_dtype);
        *out_dtype = res_dtype;
        return 0;
    }

fail:
    Py_XDECREF(*out_dtype);
    *out_dtype = NULL;
    return -1;
}

#undef RETRY_WITH_STRING
#undef RETRY_WITH_UNICODE

/* new reference */
NPY_NO_EXPORT PyArray_Descr *
_array_typedescr_fromstr(char *c_str)
{
    PyArray_Descr *descr = NULL;
    PyObject *stringobj = PyString_FromString(c_str);

    if (stringobj == NULL) {
        return NULL;
    }
    if (PyArray_DescrConverter(stringobj, &descr) != NPY_SUCCEED) {
        Py_DECREF(stringobj);
        return NULL;
    }
    Py_DECREF(stringobj);
    return descr;
}


NPY_NO_EXPORT char *
index2ptr(PyArrayObject *mp, npy_intp i)
{
    npy_intp dim0;

    if (PyArray_NDIM(mp) == 0) {
        PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
        return NULL;
    }
    dim0 = PyArray_DIMS(mp)[0];
    if (check_and_adjust_index(&i, dim0, 0, NULL) < 0)
        return NULL;
    if (i == 0) {
        return PyArray_DATA(mp);
    }
    return PyArray_BYTES(mp)+i*PyArray_STRIDES(mp)[0];
}

int NI_FindObjects(PyArrayObject* input, npy_intp max_label,
                                     npy_intp* regions)
{
    npy_intp size, jj;
    NI_Iterator ii;
    char *pi;
    NPY_BEGIN_THREADS_DEF;

    NPY_BEGIN_THREADS;

    /* get input data, size and iterator: */
    pi = (void *)PyArray_DATA(input);
    size = PyArray_SIZE(input);
    if (!NI_InitPointIterator(input, &ii))
        goto exit;
    if (PyArray_NDIM(input) > 0) {
        for (jj = 0; jj < 2 * PyArray_NDIM(input) * max_label; jj++) {
            regions[jj] = -1;
        }
    } else {
        for(jj = 0; jj < max_label; jj++)
            regions[jj] = -1;
    }
    /* iterate over all points: */
    for(jj = 0 ; jj < size; jj++) {
        switch (PyArray_TYPE(input)) {
            CASE_FIND_OBJECT_POINT(NPY_BOOL, npy_bool,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_UBYTE, npy_ubyte,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_USHORT, npy_ushort,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_UINT, npy_uint,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_ULONG, npy_ulong,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_ULONGLONG, npy_ulonglong,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_BYTE, npy_byte,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_SHORT, npy_short,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_INT, npy_int,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_LONG, npy_long,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_LONGLONG, npy_longlong,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_FLOAT, npy_float,
                                   pi, regions, input, max_label, ii);
            CASE_FIND_OBJECT_POINT(NPY_DOUBLE, npy_double,
                                   pi, regions, input, max_label, ii);
        default:
            NPY_END_THREADS;
            PyErr_SetString(PyExc_RuntimeError, "data type not supported");
            goto exit;
        }
        NI_ITERATOR_NEXT(ii, pi);
    }
 exit:
    NPY_END_THREADS;
    return PyErr_Occurred() ? 0 : 1;
}

int NI_Correlate(PyArrayObject* input, PyArrayObject* weights,
                 PyArrayObject* output, NI_ExtendMode mode,
                 double cvalue, npy_intp *origins)
{
    npy_bool *pf = NULL;
    npy_intp fsize, jj, kk, filter_size = 0, border_flag_value;
    npy_intp *offsets = NULL, *oo, size;
    NI_FilterIterator fi;
    NI_Iterator ii, io;
    char *pi, *po;
    npy_double *pw;
    npy_double *ww = NULL;
    int err = 0;
    NPY_BEGIN_THREADS_DEF;

    /* get the the footprint: */
    fsize = PyArray_SIZE(weights);
    pw = (npy_double*)PyArray_DATA(weights);
    pf = malloc(fsize * sizeof(npy_bool));
    if (!pf) {
        PyErr_NoMemory();
        goto exit;
    }
    for(jj = 0; jj < fsize; jj++) {
        if (fabs(pw[jj]) > DBL_EPSILON) {
            pf[jj] = 1;
            ++filter_size;
        } else {
            pf[jj] = 0;
        }
    }
    /* copy the weights to contiguous memory: */
    ww = malloc(filter_size * sizeof(npy_double));
    if (!ww) {
        PyErr_NoMemory();
        goto exit;
    }
    jj = 0;
    for(kk = 0; kk < fsize; kk++) {
        if (pf[kk]) {
            ww[jj++] = pw[kk];
        }
    }
    /* initialize filter offsets: */
    if (!NI_InitFilterOffsets(input, pf, PyArray_DIMS(weights), origins,
                              mode, &offsets, &border_flag_value, NULL)) {
        goto exit;
    }
    /* initialize filter iterator: */
    if (!NI_InitFilterIterator(PyArray_NDIM(input), PyArray_DIMS(weights),
                               filter_size, PyArray_DIMS(input), origins,
                               &fi)) {
        goto exit;
    }
    /* initialize input element iterator: */
    if (!NI_InitPointIterator(input, &ii))
        goto exit;
    /* initialize output element iterator: */
    if (!NI_InitPointIterator(output, &io))
        goto exit;

    NPY_BEGIN_THREADS;
    /* get data pointers an array size: */
    pi = (void *)PyArray_DATA(input);
    po = (void *)PyArray_DATA(output);
    size = PyArray_SIZE(input);
    /* iterator over the elements: */
    oo = offsets;
    for(jj = 0; jj < size; jj++) {
        double tmp = 0.0;
        switch (PyArray_TYPE(input)) {
            CASE_CORRELATE_POINT(NPY_BOOL, npy_bool,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_UBYTE, npy_ubyte,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_USHORT, npy_ushort,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_UINT, npy_uint,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_ULONG, npy_ulong,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_ULONGLONG, npy_ulonglong,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_BYTE, npy_byte,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_SHORT, npy_short,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_INT, npy_int,
                                 pi, ww, oo, filter_size, cvalue, tmp,
                                 border_flag_value);
            CASE_CORRELATE_POINT(NPY_LONG, npy_long,
                                 pi, ww, oo, filter_size, cvalue, tmp,
//.........这里部分代码省略.........

static PyObject *inverse_division_inplace(PyObject *self, PyObject *args)
{
  PyObject* a = NULL;
  PyObject* b = NULL;
  npy_intp sz = 0, i;
  npy_complex64* tmp_sp = NULL;
  npy_float32 tmp2_sp;
  npy_complex64* a_data_sp = NULL;
  npy_float32* b_data_sp = NULL;
  npy_complex128* tmp_dp = NULL;
  npy_float64 tmp2_dp;
  npy_complex128* a_data_dp = NULL;
  npy_float64* b_data_dp = NULL;
  if (!PyArg_ParseTuple(args, "OO", &a, &b))
    return NULL;
  if (!(PyArray_Check(a) && PyArray_Check(b)))
    {
      PyErr_SetString(PyExc_TypeError,"arguments must be array objects");
      return NULL;
    }
  sz = PyArray_SIZE(a);

  if (sz != PyArray_SIZE(b))
    {
      PyErr_SetString(PyExc_TypeError,"argument sizes must be equal");
      return NULL;
    }
  if ((PyArray_TYPE(a) == PyArray_COMPLEX64) && (PyArray_TYPE(b) == PyArray_FLOAT32))
    {
      a_data_sp = (npy_complex64*)PyArray_DATA(a);
      b_data_sp = (npy_float32*)PyArray_DATA(b);
      for (i=0; i<sz; ++i)
	{
	  tmp_sp = a_data_sp + i;
	  if (tmp_sp->real==0.0 || (b_data_sp[i]==0.0))
	    {
	      tmp_sp->real = tmp_sp->imag = 0.0;
	    }
	  else
	    {
	      tmp2_sp = b_data_sp[i] / (tmp_sp->real * tmp_sp->real + tmp_sp->imag * tmp_sp->imag);
	      tmp_sp->real *= tmp2_sp;
	      tmp_sp->imag *= -tmp2_sp;
	    }
	}
    }
  else if ((PyArray_TYPE(a) == PyArray_COMPLEX128) && (PyArray_TYPE(b) == PyArray_FLOAT64))
    {
      a_data_dp = (npy_complex128*)PyArray_DATA(a);
      b_data_dp = (npy_float64*)PyArray_DATA(b);
      for (i=0; i<sz; ++i)
	{
	  tmp_dp = a_data_dp + i;
	  if (tmp_dp->real==0.0 || (b_data_dp[i]==0.0))
	    {
	      tmp_dp->real = tmp_dp->imag = 0.0;
	    }
	  else
	    {
	      tmp2_dp = b_data_dp[i] / (tmp_dp->real * tmp_dp->real + tmp_dp->imag * tmp_dp->imag);
	      tmp_dp->real *= tmp2_dp;
	      tmp_dp->imag *= -tmp2_dp;
	    }
	}
    }
  else
    {
      PyErr_SetString(PyExc_TypeError,"argument types must be complex64 and float32");
      return NULL;
    }
  return Py_BuildValue("");
}

static PyObject*
libtfr_stft(PyObject *self, PyObject *args)
{
        /* arguments */
        PyObject *o = NULL;
        PyArrayObject *signal = NULL;
        PyArrayObject *signal_cast = NULL;
        PyObject *o2 = NULL;
        PyArrayObject *window = NULL;
        int step;
        int N = 0;
        int Npoints, Ntapers;

        /* output data */
        npy_intp out_shape[3];
        PyArrayObject *outdata = NULL;
        int do_complex = 0;
        double *spec;
        complex double *zspec_tmp;
        npy_cdouble *zspec;

        /* internal stuff */
        mfft *mtmh;
        double *samples = NULL;
        double *windowp;

        /* parse arguments */
        if (!PyArg_ParseTuple(args, "OOi|ii", &o, &o2, &step, &N, &do_complex))
                return NULL;
        signal = (PyArrayObject*) PyArray_FromAny(o, NULL, 1, 1, NPY_CONTIGUOUS, NULL);
        if (signal==NULL) {
                PyErr_SetString(PyExc_TypeError, "Input signal must be an ndarray");
                return NULL;
        }
        window = (PyArrayObject*) PyArray_FromAny(o2, NULL, 1, 2, NPY_CONTIGUOUS, NULL);
        if (window==NULL) {
                PyErr_SetString(PyExc_TypeError, "Window must be a 1D or 2D ndarray");
                goto fail;
        }
        /* determine dimensions of window */
        if (PyArray_NDIM(window)==1) {
                Ntapers = 1;
                Npoints = PyArray_DIM(window, 0);
        }
        else {
                Ntapers = PyArray_DIM(window, 0);
                Npoints = PyArray_DIM(window, 1);
        }

        /* coerce data to proper type */
        samples = coerce_ndarray_double(signal, &signal_cast);
        if (samples==NULL) {
                PyErr_SetString(PyExc_TypeError, "Unable to cast signal to supported data type");
                goto fail;
        }
        /* need to copy the window function b/c mtm_destroy() will demalloc it */
        if (PyArray_TYPE(window)!=NPY_DOUBLE) {
                PyErr_SetString(PyExc_TypeError, "Window function must be double precision float");
                goto fail;
        }
        windowp = malloc(Npoints * Ntapers * sizeof(double));
        memcpy(windowp, PyArray_DATA(window), Npoints * Ntapers * sizeof(double));

        /* allocate outputs and do the transform */
        if (N < 1)
                N = Npoints;
        mtmh = mtm_init(N, Npoints, Ntapers, windowp, NULL);
        out_shape[0] = N/2+1;
        if (do_complex) {
                out_shape[1] = Ntapers;
                out_shape[2] = SPEC_NFRAMES(mtmh, PyArray_SIZE(signal), step);
                outdata  = (PyArrayObject*) PyArray_ZEROS(3,out_shape,NPY_CDOUBLE,1); // fortran-order
                zspec = (npy_cdouble*) PyArray_DATA(outdata);
                //printf("output dimensions: %d, %d, %d\n", out_shape[0], out_shape[1], out_shape[2]);
                zspec_tmp = (complex double*)malloc(PyArray_SIZE(outdata) * sizeof(complex double));
                mtm_zspec(mtmh, zspec_tmp, samples, PyArray_SIZE(signal), step);
                cmplx_c99tonpy(zspec_tmp, zspec, PyArray_SIZE(outdata));
                free(zspec_tmp);
        }
        else {
                out_shape[1] = SPEC_NFRAMES(mtmh, PyArray_SIZE(signal), step);
                outdata  = (PyArrayObject*) PyArray_ZEROS(2,out_shape,NPY_DOUBLE,1); // fortran-order
                spec = (double*) PyArray_DATA(outdata);
                mtm_spec(mtmh, spec, samples, PyArray_SIZE(signal), step, 0);
        }
        mtm_destroy(mtmh);

        Py_DECREF(signal);
        Py_DECREF(window);
        Py_XDECREF(signal_cast);
        return PyArray_Return(outdata);
fail:
        Py_XDECREF(signal_cast);
        Py_XDECREF(signal);
        Py_XDECREF(window);
        Py_XDECREF(outdata);
        return NULL;
}

static PyObject * oscillAnglesOfHKLs(PyObject * self, PyObject * args)
{
  PyArrayObject *hkls, *rMat_c, *bMat,
                *vInv_s, *beamVec, *etaVec;
  PyFloatObject *chi, *wavelength;
  PyArrayObject *oangs0, *oangs1;

  int dhkls, drc, dbm, dvi, dbv, dev;
  npy_intp npts, dims[2];

  double *hkls_Ptr, chi_d,
         *rMat_c_Ptr, *bMat_Ptr, wavelen_d,
         *vInv_s_Ptr, *beamVec_Ptr, *etaVec_Ptr;
  double *oangs0_Ptr, *oangs1_Ptr;

  /* Parse arguments */
  if ( !PyArg_ParseTuple(args,"OOOOOOOO",
			 &hkls, &chi,
			 &rMat_c, &bMat, &wavelength,
			 &vInv_s, &beamVec, &etaVec)) return(NULL);
  if ( hkls    == NULL || chi == NULL ||
       rMat_c  == NULL || bMat == NULL || wavelength == NULL ||
       vInv_s  == NULL || beamVec == NULL || etaVec == NULL ) return(NULL);

  /* Verify shape of input arrays */
  dhkls = PyArray_NDIM(hkls);
  drc   = PyArray_NDIM(rMat_c);
  dbm   = PyArray_NDIM(bMat);
  dvi   = PyArray_NDIM(vInv_s);
  dbv   = PyArray_NDIM(beamVec);
  dev   = PyArray_NDIM(etaVec);
  assert( dhkls == 2 && drc == 2 && dbm == 2 &&
	  dvi   == 1 && dbv == 1 && dev == 1);

  /* Verify dimensions of input arrays */
  npts = PyArray_DIMS(hkls)[0];

  assert( PyArray_DIMS(hkls)[1]    == 3 );
  assert( PyArray_DIMS(rMat_c)[0]  == 3 && PyArray_DIMS(rMat_c)[1] == 3 );
  assert( PyArray_DIMS(bMat)[0]    == 3 && PyArray_DIMS(bMat)[1]   == 3 );
  assert( PyArray_DIMS(vInv_s)[0]  == 6 );
  assert( PyArray_DIMS(beamVec)[0] == 3 );
  assert( PyArray_DIMS(etaVec)[0]  == 3 );

  /* Allocate arrays for return values */
  dims[0] = npts; dims[1] = 3;
  oangs0 = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0);
  oangs1 = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0);

  /* Grab data pointers into various arrays */
  hkls_Ptr    = (double*)PyArray_DATA(hkls);

  chi_d       = PyFloat_AsDouble((PyObject*)chi);
  wavelen_d   = PyFloat_AsDouble((PyObject*)wavelength);

  rMat_c_Ptr  = (double*)PyArray_DATA(rMat_c);
  bMat_Ptr    = (double*)PyArray_DATA(bMat);

  vInv_s_Ptr  = (double*)PyArray_DATA(vInv_s);

  beamVec_Ptr = (double*)PyArray_DATA(beamVec);
  etaVec_Ptr  = (double*)PyArray_DATA(etaVec);

  oangs0_Ptr  = (double*)PyArray_DATA(oangs0);
  oangs1_Ptr  = (double*)PyArray_DATA(oangs1);

  /* Call the computational routine */
  oscillAnglesOfHKLs_cfunc(npts, hkls_Ptr, chi_d,
			   rMat_c_Ptr, bMat_Ptr, wavelen_d,
			   vInv_s_Ptr, beamVec_Ptr, etaVec_Ptr,
			   oangs0_Ptr, oangs1_Ptr);

  // printf("chi = %g, wavelength = %g\n",PyFloat_AsDouble((PyObject*)chi),PyFloat_AsDouble((PyObject*)wavelength));
  /*
np.ascontiguousarray(hkls),chi,rMat_c,bMat,wavelength,
                                               beamVec.flatten(),etaVec.flatten()
  */
  /* Build and return the list data structure */
  return(Py_BuildValue("OO",oangs0,oangs1));
}

/*
    Takes a list of unit reciprocal lattice vectors in crystal frame to the
    specified detector-relative frame, subject to the conditions:

    1) the reciprocal lattice vector must be able to satisfy a bragg condition
    2) the associated diffracted beam must intersect the detector plane

    Required Arguments:
    gVec_c -- (n, 3) ndarray of n reciprocal lattice vectors in the CRYSTAL FRAME
    rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components to LAB FRAME
    rMat_s -- (3, 3) ndarray, the COB taking SAMPLE FRAME components to LAB FRAME
    rMat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components to SAMPLE FRAME
    tVec_d -- (3, 1) ndarray, the translation vector connecting LAB to DETECTOR
    tVec_s -- (3, 1) ndarray, the translation vector connecting LAB to SAMPLE
    tVec_c -- (3, 1) ndarray, the translation vector connecting SAMPLE to CRYSTAL

    Outputs:
    (m, 2) ndarray containing the intersections of m <= n diffracted beams
    associated with gVecs
*/
static PyObject * gvecToDetectorXY(PyObject * self, PyObject * args)
{
  PyArrayObject *gVec_c,
                *rMat_d, *rMat_s, *rMat_c,
                *tVec_d, *tVec_s, *tVec_c, 
                *beamVec;
  PyArrayObject *result;

  int dgc, drd, drs, drc, dtd, dts, dtc, dbv;
  npy_intp npts, dims[2];

  double *gVec_c_Ptr,
         *rMat_d_Ptr, *rMat_s_Ptr, *rMat_c_Ptr,
         *tVec_d_Ptr, *tVec_s_Ptr, *tVec_c_Ptr,
         *beamVec_Ptr;
  double *result_Ptr;

  /* Parse arguments */
  if ( !PyArg_ParseTuple(args,"OOOOOOOO",
			 &gVec_c,
			 &rMat_d, &rMat_s, &rMat_c,
			 &tVec_d, &tVec_s, &tVec_c,
			 &beamVec)) return(NULL);
  if ( gVec_c  == NULL ||
       rMat_d  == NULL || rMat_s == NULL || rMat_c == NULL ||
       tVec_d  == NULL || tVec_s == NULL || tVec_c == NULL ||
       beamVec == NULL ) return(NULL);

  /* Verify shape of input arrays */
  dgc = PyArray_NDIM(gVec_c);
  drd = PyArray_NDIM(rMat_d);
  drs = PyArray_NDIM(rMat_s);
  drc = PyArray_NDIM(rMat_c);
  dtd = PyArray_NDIM(tVec_d);
  dts = PyArray_NDIM(tVec_s);
  dtc = PyArray_NDIM(tVec_c);
  dbv = PyArray_NDIM(beamVec);
  assert( dgc == 2 );
  assert( drd == 2 && drs == 2 && drc == 2 );
  assert( dtd == 1 && dts == 1 && dtc == 1 );
  assert( dbv == 1 );

  /* Verify dimensions of input arrays */
  npts = PyArray_DIMS(gVec_c)[0];

  assert( PyArray_DIMS(gVec_c)[1]  == 3 );
  assert( PyArray_DIMS(rMat_d)[0]  == 3 && PyArray_DIMS(rMat_d)[1] == 3 );
  assert( PyArray_DIMS(rMat_s)[0]  == 3 && PyArray_DIMS(rMat_s)[1] == 3 );
  assert( PyArray_DIMS(rMat_c)[0]  == 3 && PyArray_DIMS(rMat_c)[1] == 3 );
  assert( PyArray_DIMS(tVec_d)[0]  == 3 );
  assert( PyArray_DIMS(tVec_s)[0]  == 3 );
  assert( PyArray_DIMS(tVec_c)[0]  == 3 );
  assert( PyArray_DIMS(beamVec)[0] == 3 );

  /* Allocate C-style array for return data */
  // result_Ptr  = malloc(2*npts*sizeof(double));
  dims[0] = npts; dims[1] = 2;
  result = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0);

  /* Grab data pointers into various arrays */
  gVec_c_Ptr  = (double*)PyArray_DATA(gVec_c);

  rMat_d_Ptr  = (double*)PyArray_DATA(rMat_d);
  rMat_s_Ptr  = (double*)PyArray_DATA(rMat_s);
  rMat_c_Ptr  = (double*)PyArray_DATA(rMat_c);

  tVec_d_Ptr  = (double*)PyArray_DATA(tVec_d);
  tVec_s_Ptr  = (double*)PyArray_DATA(tVec_s);
  tVec_c_Ptr  = (double*)PyArray_DATA(tVec_c);

  beamVec_Ptr = (double*)PyArray_DATA(beamVec);

  result_Ptr     = (double*)PyArray_DATA(result);

  /* Call the computational routine */
  gvecToDetectorXY_cfunc(npts, gVec_c_Ptr,
			 rMat_d_Ptr, rMat_s_Ptr, rMat_c_Ptr,
			 tVec_d_Ptr, tVec_s_Ptr, tVec_c_Ptr,
			 beamVec_Ptr,
			 result_Ptr);
//.........这里部分代码省略.........

 BaseType* data() {
     return reinterpret_cast<BaseType*>PyArray_DATA(this->array_);
 }

/*
    Takes a list cartesian (x, y) pairs in the detector coordinates and calculates
    the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) pairs
    with respect to the specified beam and azimth (eta) reference directions

    Required Arguments:
    xy_det -- (n, 2) ndarray or list-like input of n detector (x, y) points
    rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components to LAB FRAME
    rMat_s -- (3, 3) ndarray, the COB taking SAMPLE FRAME components to LAB FRAME
    tVec_d -- (3, 1) ndarray, the translation vector connecting LAB to DETECTOR
    tVec_s -- (3, 1) ndarray, the translation vector connecting LAB to SAMPLE
    tVec_c -- (3, 1) ndarray, the translation vector connecting SAMPLE to CRYSTAL

    Optional Keyword Arguments:
    beamVec -- (1, 3) mdarray containing the incident beam direction components in the LAB FRAME
    etaVec  -- (1, 3) mdarray containing the reference azimuth direction components in the LAB FRAME

    Outputs:
    (n, 2) ndarray containing the (tTh, eta) pairs associated with each (x, y)
    (n, 3) ndarray containing the associated G vector directions in the LAB FRAME
    associated with gVecs
*/
static PyObject * detectorXYToGvec(PyObject * self, PyObject * args)
{
  PyArrayObject *xy_det, *rMat_d, *rMat_s,
		*tVec_d, *tVec_s, *tVec_c,
                *beamVec, *etaVec;
  PyArrayObject *tTh, *eta, *gVec_l;

  int dxy, drd, drs, dtd, dts, dtc, dbv, dev;
  npy_intp npts, dims[2];

  double *xy_Ptr, *rMat_d_Ptr, *rMat_s_Ptr,
         *tVec_d_Ptr, *tVec_s_Ptr, *tVec_c_Ptr,
         *beamVec_Ptr, *etaVec_Ptr;
  double *tTh_Ptr, *eta_Ptr, *gVec_l_Ptr;

  /* Parse arguments */
  if ( !PyArg_ParseTuple(args,"OOOOOOOO",
			 &xy_det,
			 &rMat_d, &rMat_s,
			 &tVec_d, &tVec_s, &tVec_c,
			 &beamVec, &etaVec)) return(NULL);
  if ( xy_det  == NULL || rMat_d == NULL || rMat_s == NULL ||
       tVec_d  == NULL || tVec_s == NULL || tVec_c == NULL ||
       beamVec == NULL || etaVec == NULL ) return(NULL);

  /* Verify shape of input arrays */
  dxy = PyArray_NDIM(xy_det);
  drd = PyArray_NDIM(rMat_d);
  drs = PyArray_NDIM(rMat_s);
  dtd = PyArray_NDIM(tVec_d);
  dts = PyArray_NDIM(tVec_s);
  dtc = PyArray_NDIM(tVec_c);
  dbv = PyArray_NDIM(beamVec);
  dev = PyArray_NDIM(etaVec);
  assert( dxy == 2 && drd == 2 && drs == 2 &&
	  dtd == 1 && dts == 1 && dtc == 1 &&
	  dbv == 1 && dev == 1);

  /* Verify dimensions of input arrays */
  npts = PyArray_DIMS(xy_det)[0];

  assert( PyArray_DIMS(xy_det)[1]  == 2 );
  assert( PyArray_DIMS(rMat_d)[0]  == 3 && PyArray_DIMS(rMat_d)[1] == 3 );
  assert( PyArray_DIMS(rMat_s)[0]  == 3 && PyArray_DIMS(rMat_s)[1] == 3 );
  assert( PyArray_DIMS(tVec_d)[0]  == 3 );
  assert( PyArray_DIMS(tVec_s)[0]  == 3 );
  assert( PyArray_DIMS(tVec_c)[0]  == 3 );
  assert( PyArray_DIMS(beamVec)[0] == 3 );
  assert( PyArray_DIMS(etaVec)[0]  == 3 );

  /* Allocate arrays for return values */
  dims[0] = npts; dims[1] = 3;
  gVec_l = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0);

  tTh    = (PyArrayObject*)PyArray_EMPTY(1,&npts,NPY_DOUBLE,0);
  eta    = (PyArrayObject*)PyArray_EMPTY(1,&npts,NPY_DOUBLE,0);

  /* Grab data pointers into various arrays */
  xy_Ptr      = (double*)PyArray_DATA(xy_det);
  gVec_l_Ptr  = (double*)PyArray_DATA(gVec_l);

  tTh_Ptr     = (double*)PyArray_DATA(tTh);
  eta_Ptr     = (double*)PyArray_DATA(eta);

  rMat_d_Ptr  = (double*)PyArray_DATA(rMat_d);
  rMat_s_Ptr  = (double*)PyArray_DATA(rMat_s);

  tVec_d_Ptr  = (double*)PyArray_DATA(tVec_d);
  tVec_s_Ptr  = (double*)PyArray_DATA(tVec_s);
  tVec_c_Ptr  = (double*)PyArray_DATA(tVec_c);

  beamVec_Ptr = (double*)PyArray_DATA(beamVec);
  etaVec_Ptr  = (double*)PyArray_DATA(etaVec);

  /* Call the computational routine */
  detectorXYToGvec_cfunc(npts, xy_Ptr,
			 rMat_d_Ptr, rMat_s_Ptr,
			 tVec_d_Ptr, tVec_s_Ptr, tVec_c_Ptr,
//.........这里部分代码省略.........

 void* raw_data() const { return PyArray_DATA(array_); }

PyObject * vdw(PyObject* self, PyObject *args)
{
  PyArrayObject* n_obj;
  PyArrayObject* q0_obj;
  PyArrayObject* R_obj;
  PyArrayObject* cell_obj;
  PyArrayObject* pbc_obj;
  PyArrayObject* repeat_obj;
  PyArrayObject* phi_obj;
  double ddelta;
  double dD;
  int iA;
  int iB;
  PyArrayObject* rhistogram_obj;
  double drhist;
  PyArrayObject* Dhistogram_obj;
  double dDhist;
  if (!PyArg_ParseTuple(args, "OOOOOOOddiiOdOd", &n_obj, &q0_obj, &R_obj,
                        &cell_obj, &pbc_obj, &repeat_obj,
                        &phi_obj, &ddelta, &dD, &iA, &iB,
                        &rhistogram_obj, &drhist,
                        &Dhistogram_obj, &dDhist))
    return NULL;

  int ndelta = PyArray_DIMS(phi_obj)[0];
  int nD = PyArray_DIMS(phi_obj)[1];
  const double* n = (const double*)DOUBLEP(n_obj);
  const int ni = PyArray_SIZE(n_obj);
  const double* q0 = (const double*)DOUBLEP(q0_obj);
  const double (*R)[3] = (const double (*)[3])DOUBLEP(R_obj);
  const double* cell = (const double*)DOUBLEP(cell_obj);
  const char* pbc = (const char*)(PyArray_DATA(pbc_obj));
  const long* repeat = (const long*)(PyArray_DATA(repeat_obj));
  const double (*phi)[nD] = (const double (*)[nD])DOUBLEP(phi_obj);
  double* rhistogram = (double*)DOUBLEP(rhistogram_obj);
  double* Dhistogram = (double*)DOUBLEP(Dhistogram_obj);

  int nbinsr = PyArray_DIMS(rhistogram_obj)[0];
  int nbinsD = PyArray_DIMS(Dhistogram_obj)[0];

  double energy = 0.0;
  if (repeat[0] == 0 && repeat[1] == 0 && repeat[2] == 0)
    for (int i1 = iA; i1 < iB; i1++)
      {
        const double* R1 = R[i1];
        double q01 = q0[i1];
        for (int i2 = 0; i2 <= i1; i2++)
          {
            double rr = 0.0;
            for (int c = 0; c < 3; c++)
              {
                double f = R[i2][c] - R1[c];
                if (pbc[c])
                  f = fmod(f + 1.5 * cell[c], cell[c]) - 0.5 * cell[c];
                rr += f * f;
              }
            double r = sqrt(rr);
            double d1 = r * q01;
            double d2 = r * q0[i2];
            double D = 0.5 * (d1 + d2);
            double e12 = (vdwkernel(D, d1, d2, nD, ndelta, dD, ddelta, phi) *
                          n[i1] * n[i2]);
            if (i1 == i2)
              e12 /= 2.0;
            int bin = (int)(r / drhist);
            if (bin < nbinsr)
              rhistogram[bin] += e12; 
            bin = (int)(D / dDhist);
            if (bin < nbinsD)
              Dhistogram[bin] += e12; 
            energy += e12;
          }
      }
  else
    for (int i1 = iA; i1 < iB; i1++)
      {
        const double* R1 = R[i1];
        double q01 = q0[i1];
        for (int a1 = -repeat[0]; a1 <= repeat[0]; a1++)
          for (int a2 = -repeat[1]; a2 <= repeat[1]; a2++)
            for (int a3 = -repeat[2]; a3 <= repeat[2]; a3++)
              {
                double x = 0.5;
                int i2max = ni-1;
                if (a1 == 0 && a2 == 0 && a3 == 0)
                  {
                    i2max = i1;
                    x = 1.0;
                  }
                double R1a[3] = {R1[0] + a1 * cell[0],
                                 R1[1] + a2 * cell[1],
                                 R1[2] + a3 * cell[2]};
                for (int i2 = 0; i2 <= i2max; i2++)
                  {
                    double rr = 0.0;
                    for (int c = 0; c < 3; c++)
                      {
                        double f = R[i2][c] - R1a[c];
                        rr += f * f;
                      }
//.........这里部分代码省略.........

PyObject *embedBoundsMatrix(python::object boundsMatArg, int maxIters = 10,
                            bool randomizeOnFailure = false,
                            int numZeroFail = 2,
                            python::list weights = python::list(),
                            int randomSeed = -1) {
  PyObject *boundsMatObj = boundsMatArg.ptr();
  if (!PyArray_Check(boundsMatObj))
    throw_value_error("Argument isn't an array");

  PyArrayObject *boundsMat = reinterpret_cast<PyArrayObject *>(boundsMatObj);
  // get the dimensions of the array
  unsigned int nrows = PyArray_DIM(boundsMat, 0);
  unsigned int ncols = PyArray_DIM(boundsMat, 1);
  if (nrows != ncols) throw_value_error("The array has to be square");
  if (nrows <= 0) throw_value_error("The array has to have a nonzero size");
  if (PyArray_DESCR(boundsMat)->type_num != NPY_DOUBLE)
    throw_value_error("Only double arrays are currently supported");

  unsigned int dSize = nrows * nrows;
  auto *cData = new double[dSize];
  double *inData = reinterpret_cast<double *>(PyArray_DATA(boundsMat));
  memcpy(static_cast<void *>(cData), static_cast<const void *>(inData),
         dSize * sizeof(double));

  DistGeom::BoundsMatrix::DATA_SPTR sdata(cData);
  DistGeom::BoundsMatrix bm(nrows, sdata);

  auto *positions = new RDGeom::Point3D[nrows];
  std::vector<RDGeom::Point *> posPtrs;
  for (unsigned int i = 0; i < nrows; i++) {
    posPtrs.push_back(&positions[i]);
  }

  RDNumeric::DoubleSymmMatrix distMat(nrows, 0.0);

  // ---- ---- ---- ---- ---- ---- ---- ---- ----
  // start the embedding:
  bool gotCoords = false;
  for (int iter = 0; iter < maxIters && !gotCoords; iter++) {
    // pick a random distance matrix
    DistGeom::pickRandomDistMat(bm, distMat, randomSeed);

    // and embed it:
    gotCoords = DistGeom::computeInitialCoords(
        distMat, posPtrs, randomizeOnFailure, numZeroFail, randomSeed);

    // update the seed:
    if (randomSeed >= 0) randomSeed += iter * 999;
  }

  if (gotCoords) {
    std::map<std::pair<int, int>, double> weightMap;
    unsigned int nElems = PySequence_Size(weights.ptr());
    for (unsigned int entryIdx = 0; entryIdx < nElems; entryIdx++) {
      PyObject *entry = PySequence_GetItem(weights.ptr(), entryIdx);
      if (!PySequence_Check(entry) || PySequence_Size(entry) != 3) {
        throw_value_error("weights argument must be a sequence of 3-sequences");
      }
      int idx1 = PyInt_AsLong(PySequence_GetItem(entry, 0));
      int idx2 = PyInt_AsLong(PySequence_GetItem(entry, 1));
      double w = PyFloat_AsDouble(PySequence_GetItem(entry, 2));
      weightMap[std::make_pair(idx1, idx2)] = w;
    }
    DistGeom::VECT_CHIRALSET csets;
    ForceFields::ForceField *field =
        DistGeom::constructForceField(bm, posPtrs, csets, 0.0, 0.0, &weightMap);
    CHECK_INVARIANT(field, "could not build dgeom force field");
    field->initialize();
    if (field->calcEnergy() > 1e-5) {
      int needMore = 1;
      while (needMore) {
        needMore = field->minimize();
      }
    }
    delete field;
  } else {
    throw_value_error("could not embed matrix");
  }

  // ---- ---- ---- ---- ---- ---- ---- ---- ----
  // construct the results matrix:
  npy_intp dims[2];
  dims[0] = nrows;
  dims[1] = 3;
  PyArrayObject *res = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_DOUBLE);
  double *resData = reinterpret_cast<double *>(PyArray_DATA(res));
  for (unsigned int i = 0; i < nrows; i++) {
    unsigned int iTab = i * 3;
    for (unsigned int j = 0; j < 3; ++j) {
      resData[iTab + j] = positions[i][j];  //.x;
    }
  }
  delete[] positions;

  return PyArray_Return(res);
}

static PyObject* find_changes(PyObject* self, PyObject* args) {
    PyObject* arg1 = NULL;
    PyObject* arg2 = NULL;
    if (!PyArg_ParseTuple(args, "OO", &arg1, &arg2)) return NULL;
    PyObject* np_data = PyArray_FROM_OTF(arg1, NPY_FLOAT32, NPY_IN_ARRAY);
    PyObject* np_penalties = PyArray_FROM_OTF(arg2, NPY_FLOAT32, NPY_IN_ARRAY);
    if (np_data == NULL || np_penalties == NULL) {
        Py_XDECREF(np_data);
        Py_XDECREF(np_penalties);
        return NULL;
    }
    int T = PyArray_DIM(np_data, 0);
    if (T != PyArray_DIM(np_penalties, 0) + 1) {
        PyErr_SetString(PyExc_ValueError, "Dimensions of data and penalty are not compatible");
        Py_DECREF(np_data);
        Py_DECREF(np_penalties);
        return NULL;
    }
    float* data = (float*) PyArray_DATA(np_data);
    float* penalties = (float*) PyArray_DATA(np_penalties);
    double* vals = new double[T];
    int* prev = new int[T];
    std::list<SuffStat> checks;
    checks.push_back(SuffStat(0));
    for (int t = 0; t < MIN_SEP - 1; ++t)
        checks.back().add(data[t]);
    for (int t = MIN_SEP - 1; t < T; ++t) {
        double max_val = -std::numeric_limits<double>::max();
        int max_ind = -1;
        for (std::list<SuffStat>::iterator iter = checks.begin(); iter != checks.end(); ++iter) {
            if (iter->prune_t == t) {
                iter = checks.erase(iter);
                --iter;
                continue;
            }
            iter->add(data[t]);
            double val = iter->ll();
            if (iter->t > 0)
                val += vals[iter->t - 1] - penalties[iter->t - 1];
            iter->cost = val;
            if (val > max_val) {
                max_val = val;
                max_ind = iter->t;
            }
        }
        vals[t] = max_val;
        prev[t] = max_ind;
        for (std::list<SuffStat>::iterator iter = checks.begin(); iter != checks.end(); ++iter) {
            if (t < T-1 && iter->prune_t == -1 && iter->cost < vals[t] - penalties[t])
                iter->prune_t = t + MIN_SEP;
        }
        if (t - MIN_SEP + 2 >= MIN_SEP) {
            checks.push_back(SuffStat(t-MIN_SEP+2));
            for (int s = t - MIN_SEP + 2; s <= t; ++s)
                checks.back().add(data[s]);
        }
    }
    PyObject* changes = PyList_New(0);
    int ind = prev[T-1];
    while (ind > 1) {
        PyObject* num = PyInt_FromLong(ind);
        PyList_Append(changes, num);
        Py_DECREF(num);
        ind = prev[ind-1];
    }
    PyList_Reverse(changes);
    delete[] vals;
    delete[] prev;
    Py_DECREF(np_data);
    Py_DECREF(np_penalties);
    return changes;
}

static PyObject* w_parstack(PyObject *module, PyObject *args) {

    PyObject *arrays, *offsets, *shifts, *weights, *arr;
    PyObject *result;
    int method, nparallel;
    size_t narrays, nshifts, nweights;
    size_t *clengths;
    size_t lengthout;
    int32_t offsetout;
    int lengthout_arg;
    int32_t *coffsets, *cshifts;
    double *cweights, *cresult;
    double **carrays;
    npy_intp array_dims[1];
    size_t i;
    int err;

    carrays = NULL;
    clengths = NULL;
    struct module_state *st = GETSTATE(module);

    if (!PyArg_ParseTuple(args, "OOOOiiiOi", &arrays, &offsets, &shifts,
                          &weights, &method, &lengthout_arg, &offsetout, &result, &nparallel)) {

        PyErr_SetString(
            st->error,
            "usage parstack(arrays, offsets, shifts, weights, method, lengthout, offsetout, result, nparallel)" );

        return NULL;
    }
    if (!good_array(offsets, NPY_INT32)) return NULL;
    if (!good_array(shifts, NPY_INT32)) return NULL;
    if (!good_array(weights, NPY_DOUBLE)) return NULL;
    if (result != Py_None && !good_array(result, NPY_DOUBLE)) return NULL;

    coffsets = PyArray_DATA((PyArrayObject*)offsets);
    narrays = PyArray_SIZE((PyArrayObject*)offsets);

    cshifts = PyArray_DATA((PyArrayObject*)shifts);
    nshifts = PyArray_SIZE((PyArrayObject*)shifts);

    cweights = PyArray_DATA((PyArrayObject*)weights);
    nweights = PyArray_SIZE((PyArrayObject*)weights);

    nshifts /= narrays;
    nweights /= narrays;

    if (nshifts != nweights) {
        PyErr_SetString(st->error, "weights.size != shifts.size" );
        return NULL;
    }

    if (!PyList_Check(arrays)) {
        PyErr_SetString(st->error, "arg #1 must be a list of NumPy arrays.");
        return NULL;
    }

    if ((size_t)PyList_Size(arrays) != narrays) {
        PyErr_SetString(st->error, "len(offsets) != len(arrays)");
        return NULL;
    }

    carrays = (double**)calloc(narrays, sizeof(double*));
    if (carrays == NULL) {
        PyErr_SetString(st->error, "alloc failed");
        return NULL;
    }

    clengths = (size_t*)calloc(narrays, sizeof(size_t));
    if (clengths == NULL) {
        PyErr_SetString(st->error, "alloc failed");
        free(carrays);
        return NULL;
    }

    for (i=0; i<narrays; i++) {
        arr = PyList_GetItem(arrays, i);
        if (!good_array(arr, NPY_DOUBLE)) {
            free(carrays);
            free(clengths);
            return NULL;
        }
        carrays[i] = PyArray_DATA((PyArrayObject*)arr);
        clengths[i] = PyArray_SIZE((PyArrayObject*)arr);
    }
    if (lengthout_arg < 0) {
        err = parstack_config(narrays, coffsets, clengths, nshifts, cshifts,
                              cweights, method, &lengthout, &offsetout);

        if (err != 0) {
            PyErr_SetString(st->error, "parstack_config() failed");
            free(carrays);
            free(clengths);
            return NULL;
        }
    } else {
        lengthout = (size_t)lengthout_arg;
    }

    if (method == 0) {
//.........这里部分代码省略.........

static int
array_getbuffer(PyObject *obj, Py_buffer *view, int flags)
{
    PyArrayObject *self;
    _buffer_info_t *info = NULL;

    self = (PyArrayObject*)obj;

    /* Check whether we can provide the wanted properties */
    if ((flags & PyBUF_C_CONTIGUOUS) == PyBUF_C_CONTIGUOUS &&
            !PyArray_CHKFLAGS(self, NPY_ARRAY_C_CONTIGUOUS)) {
        PyErr_SetString(PyExc_ValueError, "ndarray is not C-contiguous");
        goto fail;
    }
    if ((flags & PyBUF_F_CONTIGUOUS) == PyBUF_F_CONTIGUOUS &&
            !PyArray_CHKFLAGS(self, NPY_ARRAY_F_CONTIGUOUS)) {
        PyErr_SetString(PyExc_ValueError, "ndarray is not Fortran contiguous");
        goto fail;
    }
    if ((flags & PyBUF_ANY_CONTIGUOUS) == PyBUF_ANY_CONTIGUOUS
            && !PyArray_ISONESEGMENT(self)) {
        PyErr_SetString(PyExc_ValueError, "ndarray is not contiguous");
        goto fail;
    }
    if ((flags & PyBUF_STRIDES) != PyBUF_STRIDES &&
            !PyArray_CHKFLAGS(self, NPY_ARRAY_C_CONTIGUOUS)) {
        /* Non-strided N-dim buffers must be C-contiguous */
        PyErr_SetString(PyExc_ValueError, "ndarray is not C-contiguous");
        goto fail;
    }
    if ((flags & PyBUF_WRITEABLE) == PyBUF_WRITEABLE) {
        if (PyArray_FailUnlessWriteable(self, "buffer source array") < 0) {
            goto fail;
        }
    }
    /*
     * If a read-only buffer is requested on a read-write array, we return a
     * read-write buffer, which is dubious behavior. But that's why this call
     * is guarded by PyArray_ISWRITEABLE rather than (flags &
     * PyBUF_WRITEABLE).
     */
    if (PyArray_ISWRITEABLE(self)) {
        if (array_might_be_written(self) < 0) {
            goto fail;
        }
    }

    if (view == NULL) {
        PyErr_SetString(PyExc_ValueError, "NULL view in getbuffer");
        goto fail;
    }

    /* Fill in information */
    info = _buffer_get_info(obj);
    if (info == NULL) {
        goto fail;
    }

    view->buf = PyArray_DATA(self);
    view->suboffsets = NULL;
    view->itemsize = PyArray_ITEMSIZE(self);
    view->readonly = !PyArray_ISWRITEABLE(self);
    view->internal = NULL;
    view->len = PyArray_NBYTES(self);
    if ((flags & PyBUF_FORMAT) == PyBUF_FORMAT) {
        view->format = info->format;
    } else {
        view->format = NULL;
    }
    if ((flags & PyBUF_ND) == PyBUF_ND) {
        view->ndim = info->ndim;
        view->shape = info->shape;
    }
    else {
        view->ndim = 0;
        view->shape = NULL;
    }
    if ((flags & PyBUF_STRIDES) == PyBUF_STRIDES) {
        view->strides = info->strides;
    }
    else {
        view->strides = NULL;
    }
    view->obj = (PyObject*)self;

    Py_INCREF(self);
    return 0;

fail:
    return -1;
}