C++ PyArray_DIM函数代码示例

//.........这里部分代码省略.........
#endif
  }
  else {
    input_buffer_format = TYPE_RGBA_8;
  }

  input_buffer_to_nparray = cmsCreateTransform
        (input_buffer_profile, input_buffer_format,
         nparray_data_profile, TYPE_RGBA_8,
         INTENT_PERCEPTUAL, 0);

  width = png_get_image_width(png_ptr, info_ptr);
  height = png_get_image_height(png_ptr, info_ptr);
  rows_left = height;

  while (rows_left) {
    PyObject *pyarr = NULL;
    uint32_t rows = 0;
    uint32_t row = 0;
    const uint8_t input_buf_bytes_per_pixel = (bit_depth==8) ? 4 : 8;
    const uint32_t input_buf_row_stride = sizeof(png_byte) * width
                                          * input_buf_bytes_per_pixel;
    png_byte *input_buffer = NULL;
    png_bytep *input_buf_row_pointers = NULL;

    pyarr = PyObject_CallFunction(get_buffer_callback, "ii", width, height);
    if (! pyarr) {
      PyErr_Format(PyExc_RuntimeError, "Get-buffer callback failed");
      goto cleanup;
    }
#ifdef HEAVY_DEBUG
    //assert(PyArray_ISCARRAY(arr));
    assert(PyArray_NDIM(pyarr) == 3);
    assert(PyArray_DIM(pyarr, 1) == width);
    assert(PyArray_DIM(pyarr, 2) == 4);
    assert(PyArray_TYPE(pyarr) == NPY_UINT8);
    assert(PyArray_ISBEHAVED(pyarr));
    assert(PyArray_STRIDE(pyarr, 1) == 4*sizeof(uint8_t));
    assert(PyArray_STRIDE(pyarr, 2) ==   sizeof(uint8_t));
#endif
    rows = PyArray_DIM(pyarr, 0);

    if (rows > rows_left) {
      PyErr_Format(PyExc_RuntimeError,
                   "Attempt to read %d rows from the PNG, "
                   "but only %d are left",
                   rows, rows_left);
      goto cleanup;
    }

    input_buffer = (png_byte *) malloc(rows * input_buf_row_stride);
    input_buf_row_pointers = (png_bytep *)malloc(rows * sizeof(png_bytep));
    for (row=0; row<rows; row++) {
      input_buf_row_pointers[row] = input_buffer + (row * input_buf_row_stride);
    }

    png_read_rows(png_ptr, input_buf_row_pointers, NULL, rows);
    rows_left -= rows;

    for (row=0; row<rows; row++) {
      uint8_t *pyarr_row = (uint8_t *)PyArray_DATA(pyarr)
                         + row*PyArray_STRIDE(pyarr, 0);
      uint8_t *input_row = input_buf_row_pointers[row];
      // Really minimal fake colour management. Just remaps to sRGB.
      cmsDoTransform(input_buffer_to_nparray, input_row, pyarr_row, width);
      // lcms2 ignores alpha, so copy that verbatim

Py::Object
_path_module::point_in_path_collection(const Py::Tuple& args)
{
    args.verify_length(9);

    //segments, trans, clipbox, colors, linewidths, antialiaseds
    double                  x                = Py::Float(args[0]);
    double                  y                = Py::Float(args[1]);
    double                  radius           = Py::Float(args[2]);
    agg::trans_affine       master_transform = py_to_agg_transformation_matrix(args[3].ptr());
    Py::SeqBase<Py::Object> paths            = args[4];
    Py::SeqBase<Py::Object> transforms_obj   = args[5];
    Py::SeqBase<Py::Object> offsets_obj      = args[6];
    agg::trans_affine       offset_trans     = py_to_agg_transformation_matrix(args[7].ptr());
    bool                    filled           = Py::Int(args[8]);

    PyArrayObject* offsets = (PyArrayObject*)PyArray_FromObject(
        offsets_obj.ptr(), PyArray_DOUBLE, 0, 2);
    if (!offsets ||
            (PyArray_NDIM(offsets) == 2 && PyArray_DIM(offsets, 1) != 2) ||
            (PyArray_NDIM(offsets) == 1 && PyArray_DIM(offsets, 0) != 0))
    {
        Py_XDECREF(offsets);
        throw Py::ValueError("Offsets array must be Nx2");
    }

    size_t Npaths      = paths.length();
    size_t Noffsets    = offsets->dimensions[0];
    size_t N           = std::max(Npaths, Noffsets);
    size_t Ntransforms = std::min(transforms_obj.length(), N);
    size_t i;

    // Convert all of the transforms up front
    typedef std::vector<agg::trans_affine> transforms_t;
    transforms_t transforms;
    transforms.reserve(Ntransforms);
    for (i = 0; i < Ntransforms; ++i)
    {
        agg::trans_affine trans = py_to_agg_transformation_matrix
                                  (transforms_obj[i].ptr(), false);
        trans *= master_transform;
        transforms.push_back(trans);
    }

    Py::List result;
    agg::trans_affine trans;

    for (i = 0; i < N; ++i)
    {
        PathIterator path(paths[i % Npaths]);

        if (Ntransforms)
        {
            trans = transforms[i % Ntransforms];
        }
        else
        {
            trans = master_transform;
        }

        if (Noffsets)
        {
            double xo = *(double*)PyArray_GETPTR2(offsets, i % Noffsets, 0);
            double yo = *(double*)PyArray_GETPTR2(offsets, i % Noffsets, 1);
            offset_trans.transform(&xo, &yo);
            trans *= agg::trans_affine_translation(xo, yo);
        }

        if (filled)
        {
            if (::point_in_path(x, y, path, trans))
                result.append(Py::Int((int)i));
        }
        else
        {
            if (::point_on_path(x, y, radius, path, trans))
                result.append(Py::Int((int)i));
        }
    }

    return result;
}

Py::Object
_path_module::update_path_extents(const Py::Tuple& args)
{
    args.verify_length(5);

    double x0, y0, x1, y1;
    PathIterator path(args[0]);
    agg::trans_affine trans = py_to_agg_transformation_matrix(
        args[1].ptr(), false);

    if (!py_convert_bbox(args[2].ptr(), x0, y0, x1, y1))
    {
        throw Py::ValueError(
            "Must pass Bbox object as arg 3 of update_path_extents");
    }
    Py::Object minpos_obj = args[3];
    bool ignore = Py::Boolean(args[4]);

    double xm, ym;
    PyArrayObject* input_minpos = NULL;
    try
    {
        input_minpos = (PyArrayObject*)PyArray_FromObject(
            minpos_obj.ptr(), PyArray_DOUBLE, 1, 1);
        if (!input_minpos || PyArray_DIM(input_minpos, 0) != 2)
        {
            throw Py::TypeError(
                "Argument 4 to update_path_extents must be a length-2 numpy array.");
        }
        xm = *(double*)PyArray_GETPTR1(input_minpos, 0);
        ym = *(double*)PyArray_GETPTR1(input_minpos, 1);
    }
    catch (...)
    {
        Py_XDECREF(input_minpos);
        throw;
    }
    Py_XDECREF(input_minpos);

    npy_intp extent_dims[] = { 2, 2, 0 };
    double* extents_data = NULL;
    npy_intp minpos_dims[] = { 2, 0 };
    double* minpos_data = NULL;
    PyArrayObject* extents = NULL;
    PyArrayObject* minpos = NULL;
    bool changed = false;

    try
    {
        extents = (PyArrayObject*)PyArray_SimpleNew
                  (2, extent_dims, PyArray_DOUBLE);
        if (extents == NULL)
        {
            throw Py::MemoryError("Could not allocate result array");
        }
        minpos = (PyArrayObject*)PyArray_SimpleNew
                 (1, minpos_dims, PyArray_DOUBLE);
        if (minpos == NULL)
        {
            throw Py::MemoryError("Could not allocate result array");
        }

        extents_data = (double*)PyArray_DATA(extents);
        minpos_data = (double*)PyArray_DATA(minpos);

        if (ignore)
        {
            extents_data[0] = std::numeric_limits<double>::infinity();
            extents_data[1] = std::numeric_limits<double>::infinity();
            extents_data[2] = -std::numeric_limits<double>::infinity();
            extents_data[3] = -std::numeric_limits<double>::infinity();
            minpos_data[0] = std::numeric_limits<double>::infinity();
            minpos_data[1] = std::numeric_limits<double>::infinity();
        }
        else
        {
            if (x0 > x1)
            {
                extents_data[0] = std::numeric_limits<double>::infinity();
                extents_data[2] = -std::numeric_limits<double>::infinity();
            }
            else
            {
                extents_data[0] = x0;
                extents_data[2] = x1;
            }
            if (y0 > y1)
            {
                extents_data[1] = std::numeric_limits<double>::infinity();
                extents_data[3] = -std::numeric_limits<double>::infinity();
            }
            else
            {
                extents_data[1] = y0;
                extents_data[3] = y1;
            }
            minpos_data[0] = xm;
            minpos_data[1] = ym;
        }

//.........这里部分代码省略.........

PyObject *
py_pair_distribution(PyObject *self, PyObject *args)
{
  PyObject *i_arr, *r_arr;
  int nbins;
  double cutoff;

  if (!PyArg_ParseTuple(args, "O!O!id", &PyArray_Type, &i_arr,
			&PyArray_Type, &r_arr, &nbins, &cutoff))
    return NULL;

  if (PyArray_NDIM(i_arr) != 1 || PyArray_TYPE(i_arr) != NPY_INT) {
    PyErr_SetString(PyExc_TypeError, "First argument needs to be "
                    "one-dimensional integer array.");
    return NULL;
  }
  if (PyArray_NDIM(r_arr) != 1 || PyArray_TYPE(r_arr) != NPY_DOUBLE) {
    PyErr_SetString(PyExc_TypeError, "Second argument needs to be "
                    "one-dimensional double array.");
    return NULL;
  }

  npy_intp npairs = PyArray_DIM(i_arr, 0);
  if (PyArray_DIM(r_arr, 0) != npairs) {
    PyErr_SetString(PyExc_RuntimeError, "First two arguments need to be arrays "
                    "of identical length.");
    return NULL;
  }

  npy_intp dim = nbins;
  PyObject *h_arr = PyArray_ZEROS(1, &dim, NPY_DOUBLE, 1);
  PyObject *h2_arr = PyArray_ZEROS(1, &dim, NPY_DOUBLE, 1);
  PyObject *tmp_arr = PyArray_ZEROS(1, &dim, NPY_INT, 1);

  npy_int *i = PyArray_DATA(i_arr);
  double *r = PyArray_DATA(r_arr);
  double *h = PyArray_DATA(h_arr);
  double *h2 = PyArray_DATA(h2_arr);
  npy_int *tmp = PyArray_DATA(tmp_arr);

  npy_int last_i = i[0];
  memset(tmp, 0, nbins*sizeof(npy_int));
  int nat = 1, p;
  for (p = 0; p < npairs; p++) {
    if (last_i != i[p]) {
      int bin;
      for (bin = 0; bin < nbins; bin++) {
	h[bin] += tmp[bin];
	h2[bin] += tmp[bin]*tmp[bin];
      }
      memset(tmp, 0, nbins*sizeof(npy_int));
      last_i = i[p];
      nat++;
    }

    int bin = (int) (nbins*r[p]/cutoff);
    if (bin >= 0 && bin < nbins) {
      tmp[bin]++;
    }
  }
  int bin;
  for (bin = 0; bin < nbins; bin++) {
    h[bin] += tmp[bin];
    h2[bin] += tmp[bin]*tmp[bin];

    double r1 = bin*cutoff/nbins, r2 = (bin+1)*cutoff/nbins;
    double binvol = 4*M_PI/3*(r2*r2*r2 - r1*r1*r1);

    h[bin] /= nat*binvol;
    h2[bin] /= nat*binvol*binvol;
    h2[bin] -= h[bin]*h[bin];
  }

  Py_DECREF(tmp_arr);

  return Py_BuildValue("OO", h_arr, h2_arr);
}

Py::Object TriModule::new_triangulation(const Py::Tuple &args)
{
    _VERBOSE("TriModule::new_triangulation");
    args.verify_length(6);

    // x and y.
    PyArrayObject* x = (PyArrayObject*)PyArray_ContiguousFromObject(
                           args[0].ptr(), PyArray_DOUBLE, 1, 1);
    PyArrayObject* y = (PyArrayObject*)PyArray_ContiguousFromObject(
                           args[1].ptr(), PyArray_DOUBLE, 1, 1);
    if (x == 0 || y == 0 || PyArray_DIM(x,0) != PyArray_DIM(y,0)) {
        Py_XDECREF(x);
        Py_XDECREF(y);
        throw Py::ValueError("x and y must be 1D arrays of the same length");
    }

    // triangles.
    PyArrayObject* triangles = (PyArrayObject*)PyArray_ContiguousFromObject(
                                   args[2].ptr(), PyArray_INT, 2, 2);
    if (triangles == 0 || PyArray_DIM(triangles,1) != 3) {
        Py_XDECREF(x);
        Py_XDECREF(y);
        Py_XDECREF(triangles);
        throw Py::ValueError("triangles must be a 2D array of shape (?,3)");
    }

    // Optional mask.
    PyArrayObject* mask = 0;
    if (args[3].ptr() != 0 && args[3] != Py::None())
    {
        mask = (PyArrayObject*)PyArray_ContiguousFromObject(
                   args[3].ptr(), PyArray_BOOL, 1, 1);
        if (mask == 0 || PyArray_DIM(mask,0) != PyArray_DIM(triangles,0)) {
            Py_XDECREF(x);
            Py_XDECREF(y);
            Py_XDECREF(triangles);
            Py_XDECREF(mask);
            throw Py::ValueError(
                "mask must be a 1D array with the same length as the triangles array");
        }
    }

    // Optional edges.
    PyArrayObject* edges = 0;
    if (args[4].ptr() != 0 && args[4] != Py::None())
    {
        edges = (PyArrayObject*)PyArray_ContiguousFromObject(
                    args[4].ptr(), PyArray_INT, 2, 2);
        if (edges == 0 || PyArray_DIM(edges,1) != 2) {
            Py_XDECREF(x);
            Py_XDECREF(y);
            Py_XDECREF(triangles);
            Py_XDECREF(mask);
            Py_XDECREF(edges);
            throw Py::ValueError("edges must be a 2D array with shape (?,2)");
        }
    }

    // Optional neighbors.
    PyArrayObject* neighbors = 0;
    if (args[5].ptr() != 0 && args[5] != Py::None())
    {
        neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(
                        args[5].ptr(), PyArray_INT, 2, 2);
        if (neighbors == 0 ||
            PyArray_DIM(neighbors,0) != PyArray_DIM(triangles,0) ||
            PyArray_DIM(neighbors,1) != PyArray_DIM(triangles,1)) {
            Py_XDECREF(x);
            Py_XDECREF(y);
            Py_XDECREF(triangles);
            Py_XDECREF(mask);
            Py_XDECREF(edges);
            Py_XDECREF(neighbors);
            throw Py::ValueError(
                "neighbors must be a 2D array with the same shape as the triangles array");
        }
    }

    return Py::asObject(new Triangulation(x, y, triangles, mask, edges, neighbors));
}

static PyObject* ccdToQ(PyObject *self, PyObject *args, PyObject *kwargs){
  static char *kwlist[] = { "angles", "mode", "ccd_size", "ccd_pixsize", 
			    "ccd_cen", "dist", "wavelength", 
			    "UBinv", "outarray", NULL };
  PyObject *angles = NULL;
  PyObject *_angles = NULL;
  PyObject *_ubinv = NULL;
  PyObject *ubinv = NULL;
  PyObject *_outarray = NULL;
  PyObject *qOut = NULL;
  CCD ccd;
  npy_intp dims[2];
  npy_intp nimages;
  int i, j, t, stride;
  int ndelgam;
  int mode;

  _float lambda;

  _float *anglesp;
  _float *qOutp;
  _float *ubinvp;
  _float UBI[3][3];

#ifdef USE_THREADS
  pthread_t thread[NTHREADS];
  int iret[NTHREADS];
#endif
  imageThreadData threadData[NTHREADS];

  if(!PyArg_ParseTupleAndKeywords(args, kwargs, "Oi(ii)(dd)(dd)ddO|O", kwlist,
				  &_angles,
				  &mode,
				  &ccd.xSize, &ccd.ySize,
				  &ccd.xPixSize, &ccd.yPixSize, 
				  &ccd.xCen, &ccd.yCen,
				  &ccd.dist,
				  &lambda,
				  &_ubinv,
				  &_outarray)){
    return NULL;
  }

  angles = PyArray_FROMANY(_angles, NPY_DOUBLE, 2, 2, NPY_IN_ARRAY);
  if(!angles){
    PyErr_SetString(PyExc_ValueError, "angles must be a 2-D array of floats");
    goto cleanup;
  }
  
  ubinv = PyArray_FROMANY(_ubinv, NPY_DOUBLE, 2, 2, NPY_IN_ARRAY);
  if(!ubinv){
    PyErr_SetString(PyExc_ValueError, "ubinv must be a 2-D array of floats");
    goto cleanup;
  }

  ubinvp = (_float *)PyArray_DATA(ubinv);
  for(i=0;i<3;i++){
    UBI[i][0] = -1.0 * ubinvp[2];
    UBI[i][1] = ubinvp[1];
    UBI[i][2] = ubinvp[0];
    ubinvp+=3;
  }
  
  nimages = PyArray_DIM(angles, 0);
  ndelgam = ccd.xSize * ccd.ySize;

  dims[0] = nimages * ndelgam;
  dims[1] = 4;
  if(!_outarray){
    // Create new numpy array
    // fprintf(stderr, "**** Creating new array\n");
    qOut = PyArray_SimpleNew(2, dims, NPY_DOUBLE);
    if(!qOut){
      goto cleanup;
    }
  } else {
    qOut = PyArray_FROMANY(_outarray, NPY_DOUBLE, 2, 2, NPY_INOUT_ARRAY);
    if(!qOut){
      PyErr_SetString(PyExc_ValueError, "outarray must be a 2-D array of floats");
      goto cleanup;
    }
    if(PyArray_Size(qOut) != (4 * nimages * ndelgam)){
      PyErr_SetString(PyExc_ValueError, "outarray is of the wrong size");
      goto cleanup;
    }
  }
  anglesp = (_float *)PyArray_DATA(angles);
  qOutp = (_float *)PyArray_DATA(qOut);

  stride = nimages / NTHREADS;
  for(t=0;t<NTHREADS;t++){
    // Setup threads
    // Allocate memory for delta/gamma pairs
    
    threadData[t].ccd = &ccd;
    threadData[t].anglesp = anglesp;
    threadData[t].qOutp = qOutp;
    threadData[t].ndelgam = ndelgam;
    threadData[t].lambda = lambda;
    threadData[t].mode = mode;
//.........这里部分代码省略.........

int array_levinson_nd(PyArrayObject *arr, long order,
                      PyArrayObject** alpccoeff,
                      PyArrayObject **klpccoeff, PyArrayObject **elpc)
{
	double *acoeff, *kcoeff, *tmp;
	double *err;
	double *data;
	npy_int rank;
	npy_intp alpc_size[NPY_MAXDIMS];
	npy_intp klpc_size[NPY_MAXDIMS];
	npy_intp elpc_size[NPY_MAXDIMS];
	npy_int n, nrepeat;
	int i;

	rank = PyArray_NDIM(arr);
	if (rank < 2) {
		return -1;
	}

	nrepeat = 1;
	for (i = 0; i < rank - 1; ++i) {
		nrepeat *= PyArray_DIM(arr, i);
		alpc_size[i] = PyArray_DIM(arr, i);
		klpc_size[i] = PyArray_DIM(arr, i);
		elpc_size[i] = PyArray_DIM(arr, i);
	}
	alpc_size[rank-1] = order + 1;
	klpc_size[rank-1] = order;

	*alpccoeff = (PyArrayObject*)PyArray_SimpleNew(rank, alpc_size,
                                                       PyArray_DOUBLE);
        if(*alpccoeff == NULL) {
                return -1;
        }

	*klpccoeff = (PyArrayObject*)PyArray_SimpleNew(rank, klpc_size,
						       NPY_DOUBLE);
        if(*klpccoeff == NULL) {
                goto clean_alpccoeff;
        }

        *elpc = (PyArrayObject*)PyArray_SimpleNew(rank-1, elpc_size, NPY_DOUBLE);
        if(*elpc == NULL) {
                goto clean_klpccoeff;
        }

	tmp = malloc(sizeof(*tmp) * order);
	if (tmp == NULL) {
                goto clean_elpc;
	}

	data = (double*)arr->data;
	acoeff = (double*)((*alpccoeff)->data);
	kcoeff = (double*)((*klpccoeff)->data);
	err = (double*)((*elpc)->data);
	n = PyArray_DIM(arr, rank-1);
	for(i = 0; i < nrepeat; ++i) {
		levinson(data, order, acoeff, err, kcoeff, tmp);
		data += n;
		acoeff += order + 1;
		kcoeff += order;
		err += 1;
	}

        free(tmp);
        return 0;

clean_elpc:
        Py_DECREF(*elpc);
clean_klpccoeff:
        Py_DECREF(*klpccoeff);
clean_alpccoeff:
        Py_DECREF(*alpccoeff);
	return -1;
}

static PyObject *cs_gamma_findzofA(PyObject *self, PyObject *args)
{
  PyArrayObject *Numpy_amp;
  PyObject *Numpy_zofA;
  double Gmu, alpha, *zofA, *amp;
  unsigned long int Namp;
  (void)self;	/* silence unused parameter warning */

  double z_min = 1e-20, z_max = 1e10;
  double dlnz = 0.05;
  unsigned numz = floor( (log(z_max) - log(z_min)) / dlnz );
  unsigned long int i;
  cs_cosmo_functions_t cosmofns;
  double *fz,*z;
  double a;
  gsl_interp *zofa_interp;
  gsl_interp_accel *acc_zofa = gsl_interp_accel_alloc();

  if (!PyArg_ParseTuple(args, "ddO!", &Gmu, &alpha, &PyArray_Type, &Numpy_amp))
    return NULL;

  Numpy_amp = PyArray_GETCONTIGUOUS(Numpy_amp);
  if(!Numpy_amp)
    return NULL;
  Namp = PyArray_DIM(Numpy_amp, 0);
  amp = PyArray_DATA(Numpy_amp);

  {
  npy_intp dims[1] = {Namp};
  Numpy_zofA = PyArray_SimpleNew(1, dims, NPY_DOUBLE);
  }
  zofA = PyArray_DATA((PyArrayObject *) Numpy_zofA);

  cosmofns = XLALCSCosmoFunctionsAlloc( z_min, dlnz, numz );

  zofa_interp = gsl_interp_alloc (gsl_interp_linear, cosmofns.n);

  fz = calloc( cosmofns.n, sizeof( *fz ) );
  z = calloc( cosmofns.n, sizeof( *z ) );

  /* first compute the function that relates A and z */
  /* invert order; b/c fz is a monotonically decreasing func of z */
  for ( i = cosmofns.n ; i > 0; i-- )
    {
      unsigned long int j = cosmofns.n - i;
      z[j] = cosmofns.z[i-1];
      fz[j] = pow(cosmofns.phit[i-1], 2.0/3.0) * pow(1+z[j], -1.0/3.0) / cosmofns.phiA[i-1];
    }

  gsl_interp_init (zofa_interp, fz, z, cosmofns.n);

  /* now compute the amplitudes (suitably multiplied) that are equal to fz for some z*/
  for ( i = 0; i < Namp; i++ )
    {
      a = amp[i] * pow(H0,-1.0/3.0) * pow(alpha,-2.0/3.0) / Gmu;
      /* evaluate z(fz) at fz=a */
      zofA[i] = gsl_interp_eval (zofa_interp, fz, z, a, acc_zofa );
      if(gsl_isnan(zofA[i])) {
        Py_DECREF(Numpy_zofA);
        Numpy_zofA = NULL;
        break;
      }
    }

  XLALCSCosmoFunctionsFree( cosmofns );
  Py_DECREF(Numpy_amp);
  free(fz);
  free(z);
  gsl_interp_free (zofa_interp);
  gsl_interp_accel_free(acc_zofa);

  return Numpy_zofA;
}

/*
 * digitize(x, bins, right=False) returns an array of integers the same length
 * as x. The values i returned are such that bins[i - 1] <= x < bins[i] if
 * bins is monotonically increasing, or bins[i - 1] > x >= bins[i] if bins
 * is monotonically decreasing.  Beyond the bounds of bins, returns either
 * i = 0 or i = len(bins) as appropriate. If right == True the comparison
 * is bins [i - 1] < x <= bins[i] or bins [i - 1] >= x > bins[i]
 */
NPY_NO_EXPORT PyObject *
arr_digitize(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
{
    PyObject *obj_x = NULL;
    PyObject *obj_bins = NULL;
    PyArrayObject *arr_x = NULL;
    PyArrayObject *arr_bins = NULL;
    PyObject *ret = NULL;
    npy_intp len_bins;
    int monotonic, right = 0;
    NPY_BEGIN_THREADS_DEF

    static char *kwlist[] = {"x", "bins", "right", NULL};

    if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO|i", kwlist,
                                     &obj_x, &obj_bins, &right)) {
        goto fail;
    }

    /* PyArray_SearchSorted will make `x` contiguous even if we don't */
    arr_x = (PyArrayObject *)PyArray_FROMANY(obj_x, NPY_DOUBLE, 0, 0,
                                             NPY_ARRAY_CARRAY_RO);
    if (arr_x == NULL) {
        goto fail;
    }

    /* TODO: `bins` could be strided, needs change to check_array_monotonic */
    arr_bins = (PyArrayObject *)PyArray_FROMANY(obj_bins, NPY_DOUBLE, 1, 1,
                                               NPY_ARRAY_CARRAY_RO);
    if (arr_bins == NULL) {
        goto fail;
    }

    len_bins = PyArray_SIZE(arr_bins);
    if (len_bins == 0) {
        PyErr_SetString(PyExc_ValueError, "bins must have non-zero length");
        goto fail;
    }

    NPY_BEGIN_THREADS_THRESHOLDED(len_bins)
    monotonic = check_array_monotonic((const double *)PyArray_DATA(arr_bins),
                                      len_bins);
    NPY_END_THREADS

    if (monotonic == 0) {
        PyErr_SetString(PyExc_ValueError,
                        "bins must be monotonically increasing or decreasing");
        goto fail;
    }

    /* PyArray_SearchSorted needs an increasing array */
    if (monotonic == - 1) {
        PyArrayObject *arr_tmp = NULL;
        npy_intp shape = PyArray_DIM(arr_bins, 0);
        npy_intp stride = -PyArray_STRIDE(arr_bins, 0);
        void *data = (void *)(PyArray_BYTES(arr_bins) - stride * (shape - 1));

        arr_tmp = (PyArrayObject *)PyArray_New(&PyArray_Type, 1, &shape,
                                               NPY_DOUBLE, &stride, data, 0,
                                               PyArray_FLAGS(arr_bins), NULL);
        if (!arr_tmp) {
            goto fail;
        }

        if (PyArray_SetBaseObject(arr_tmp, (PyObject *)arr_bins) < 0) {

            Py_DECREF(arr_tmp);
            goto fail;
        }
        arr_bins = arr_tmp;
    }

    ret = PyArray_SearchSorted(arr_bins, (PyObject *)arr_x,
                               right ? NPY_SEARCHLEFT : NPY_SEARCHRIGHT, NULL);
    if (!ret) {
        goto fail;
    }

    /* If bins is decreasing, ret has bins from end, not start */
    if (monotonic == -1) {
        npy_intp *ret_data =
                        (npy_intp *)PyArray_DATA((PyArrayObject *)ret);
        npy_intp len_ret = PyArray_SIZE((PyArrayObject *)ret);

        NPY_BEGIN_THREADS_THRESHOLDED(len_ret)
        while (len_ret--) {
            *ret_data = len_bins - *ret_data;
            ret_data++;
        }
        NPY_END_THREADS
    }

static PyObject *sky_map_toa_phoa_snr(
    PyObject *NPY_UNUSED(module), PyObject *args, PyObject *kwargs)
{
    /* Input arguments */
    long nside = -1;
    long npix;
    double min_distance;
    double max_distance;
    int prior_distance_power;
    double gmst;
    unsigned int nifos;
    unsigned long nsamples = 0;
    double sample_rate;
    PyObject *acors_obj;
    PyObject *responses_obj;
    PyObject *locations_obj;
    PyObject *horizons_obj;
    PyObject *toas_obj;
    PyObject *phoas_obj;
    PyObject *snrs_obj;

    /* Names of arguments */
    static const char *keywords[] = {"min_distance", "max_distance",
        "prior_distance_power", "gmst", "sample_rate", "acors", "responses",
        "locations", "horizons", "toas", "phoas", "snrs", "nside", NULL};

    /* Parse arguments */
    if (!PyArg_ParseTupleAndKeywords(args, kwargs, "ddiddOOOOOOO|l",
        keywords, &min_distance, &max_distance, &prior_distance_power, &gmst,
        &sample_rate, &acors_obj, &responses_obj, &locations_obj,
        &horizons_obj, &toas_obj, &phoas_obj, &snrs_obj, &nside)) return NULL;

    /* Determine HEALPix resolution, if specified */
    if (nside == -1)
    {
        npix = -1;
    } else {
        npix = nside2npix(nside);
        if (npix == -1)
        {
            PyErr_SetString(PyExc_ValueError, "nside must be a power of 2");
            return NULL;
        }
    }

    /* Determine number of detectors */
    {
        Py_ssize_t n = PySequence_Length(acors_obj);
        if (n < 0) return NULL;
        nifos = n;
    }

    /* Return value */
    PyObject *out = NULL;

    /* Numpy array objects */
    PyArrayObject *acors_npy[nifos], *responses_npy[nifos],
        *locations_npy[nifos], *horizons_npy = NULL, *toas_npy = NULL,
        *phoas_npy = NULL, *snrs_npy = NULL;
    memset(acors_npy, 0, sizeof(acors_npy));
    memset(responses_npy, 0, sizeof(responses_npy));
    memset(locations_npy, 0, sizeof(locations_npy));

    /* Arrays of pointers for inputs with multiple dimensions */
    const double complex *acors[nifos];
    const float (*responses[nifos])[3];
    const double *locations[nifos];

    /* Gather C-aligned arrays from Numpy types */
    INPUT_LIST_OF_ARRAYS(acors, NPY_CDOUBLE, 1,
        npy_intp dim = PyArray_DIM(npy, 0);
        if (iifo == 0)
            nsamples = dim;
        else if ((unsigned long)dim != nsamples)
        {
            PyErr_SetString(PyExc_ValueError,
                "expected elements of acors to be vectors of the same length");
            goto fail;
        }
    )

static PyObject* py_fitexpsin(PyObject *obj, PyObject *args, PyObject *kwds)
{
    PyArrayObject *data = NULL;
    PyArrayObject *fitt = NULL;
    PyArrayObject *rslt = NULL;
    PyArrayIterObject *data_it = NULL;
    PyArrayIterObject *fitt_it = NULL;
    PyArrayIterObject *rslt_it = NULL;
    Py_ssize_t newshape[NPY_MAXDIMS];
    double *poly = NULL;
    double *coef = NULL;
    double *buff = NULL;
    int i, j, error, lastaxis, numdata;
    int startcoef = -1;
    int numcoef = MAXCOEF;
    int axis = NPY_MAXDIMS;
    double deltat = 1.0;
    static char *kwlist[] = {"data", "numcoef",
                             "deltat", "axis", NULL};

    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&|idO&", kwlist,
        PyConverter_AnyDoubleArray, &data, &numcoef, &deltat,
        PyArray_AxisConverter, &axis)) return NULL;

    if (axis < 0) {
        axis += PyArray_NDIM(data);
    }
    if ((axis < 0) || (axis >= NPY_MAXDIMS)) {
        PyErr_Format(PyExc_ValueError, "invalid axis");
        goto _fail;
    }
    lastaxis = PyArray_NDIM(data) - 1;

    if ((numcoef < 1) || (numcoef > MAXCOEF)) {
        PyErr_Format(PyExc_ValueError, "numcoef out of bounds");
        goto _fail;
    }

    if (startcoef < 0) { /* start regression away from zero coefficients */
        startcoef = 4;
    }
    if (startcoef > numcoef - 2) {
        PyErr_Format(PyExc_ValueError, "startcoef out of bounds");
        goto _fail;
    }

    numdata = (int)PyArray_DIM(data, axis);
    if (numcoef > numdata)
        numcoef = numdata;

    if ((numcoef - startcoef - 1) < 3) {
        PyErr_Format(PyExc_ValueError,
            "number of coefficients insufficient to fit data");
        goto _fail;
    }

    /* fitted data */
    fitt = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(data),
                                              PyArray_DIMS(data), NPY_DOUBLE);
    if (fitt == NULL) {
        PyErr_Format(PyExc_MemoryError, "unable to allocate fitt array");
        goto _fail;
    }

    /* fitted parameters */
    j = 0;
    for (i = 0; i < PyArray_NDIM(data); i++) {
        if (i != axis)
            newshape[j++] = PyArray_DIM(data, i);
    }
    newshape[j] = 5;
    rslt = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(data),
                                              newshape, NPY_DOUBLE);
    if (rslt == NULL) {
        PyErr_Format(PyExc_MemoryError, "unable to allocate rslt array");
        goto _fail;
    }

    /* working buffer */
    buff = (double *)PyMem_Malloc(3*numdata * sizeof(double));
    if (buff == NULL) {
        PyErr_Format(PyExc_MemoryError, "unable to allocate buff array");
        goto _fail;
    }

    /* buffer for differential coefficients */
    coef = (double *)PyMem_Malloc((3+1)*(numcoef+1) * sizeof(double));
    if (coef == NULL) {
        PyErr_Format(PyExc_MemoryError, "unable to allocate coef array");
        goto _fail;
    }

    /* precalculate normalized Chebyshev polynomial */
    poly = (double *)PyMem_Malloc(numdata * (numcoef+1) * sizeof(double));
    if (poly == NULL) {
        PyErr_Format(PyExc_MemoryError, "unable to allocate poly");
        goto _fail;
    }

    error = chebypoly(numdata, numcoef, poly, 1);
//.........这里部分代码省略.........

    const double *locations[nifos];

    /* Gather C-aligned arrays from Numpy types */
    INPUT_LIST_OF_ARRAYS(acors, NPY_CDOUBLE, 1,
        npy_intp dim = PyArray_DIM(npy, 0);
        if (iifo == 0)
            nsamples = dim;
        else if ((unsigned long)dim != nsamples)
        {
            PyErr_SetString(PyExc_ValueError,
                "expected elements of acors to be vectors of the same length");
            goto fail;
        }
    )
    INPUT_LIST_OF_ARRAYS(responses, NPY_FLOAT, 2,
        if (PyArray_DIM(npy, 0) != 3 || PyArray_DIM(npy, 1) != 3)
        {
            PyErr_SetString(PyExc_ValueError,
                "expected elements of responses to be 3x3 arrays");
            goto fail;
        }
    )
    INPUT_LIST_OF_ARRAYS(locations, NPY_DOUBLE, 1,
        if (PyArray_DIM(npy, 0) != 3)
        {
            PyErr_SetString(PyExc_ValueError,
                "expected elements of locations to be vectors of length 3");
            goto fail;
        }
    )
    INPUT_VECTOR_DOUBLE_NIFOS(horizons)

static PyObject* opnorm_opnorm(PyObject *self, PyObject *args, PyObject *kwargs)
{
    PyArrayObject *A;
    double p, q, eps = 1e-10;
    size_t fifomax = 0;
    static char *kwlist[] = {"A", "p", "q", "eps", "fifomax", NULL};

    /* get arguments */

    if (! PyArg_ParseTupleAndKeywords(args, kwargs, "Odd|dn", kwlist,
                                      &A, &p, &q, &eps, &fifomax) )
    {
        PyErr_SetString(PyExc_ValueError, "error parsing arguments");
        return NULL;
    }

    /* check and retrieve dimension */

    if (PyArray_NDIM(A) != 2)
    {
        PyErr_SetString(PyExc_ValueError, "First argument must be 2D array");
        return NULL;
    }

    size_t
    m = PyArray_DIM(A, 0),
    n = PyArray_DIM(A, 1);

    /* get matrix data */

    const double *Adat = PyArray_DATA(A);

    /* create array object for maximising vector */

    PyArrayObject *vmax;
    npy_intp dims[1] = {n};

    if (! (vmax = (PyArrayObject*)PyArray_SimpleNew(1, dims, NPY_DOUBLE)))
    {
        PyErr_SetString(PyExc_RuntimeError, "failed to create output array");
        return NULL;
    }

    /* get underlying data for vmax */

    double *vdat = (double*)PyArray_DATA(vmax);

    /* run opnorm */

    double N;
    opnorm_opt_t opt = { .eps = eps, .fifomax = fifomax };
    opnorm_stats_t stats;

    int err = opnorm(Adat, row_major, m, n, p, q, opt, &N, vdat, &stats);

    if (err)
    {
        const char *msg = opnorm_strerror(err);
        switch (err)
        {
        case OPNORM_EDOM_P:
        case OPNORM_EDOM_Q:
        case OPNORM_EDOM_EPS:
            PyErr_SetString(PyExc_ValueError, msg);
            break;
        default:
            PyErr_SetString(PyExc_RuntimeError, msg);
        }
        return NULL;
    }

    return Py_BuildValue("fO{s:k,s:k,s:k,s:I}", N, vmax,
                         "neval",   stats.neval,
                         "nfifo",   stats.nfifo,
                         "fifomax", stats.fifomax,
                         "nthread", stats.nthread);
}

PyObject *
save_png_fast_progressive (char *filename,
                           int w, int h,
                           bool has_alpha,
                           PyObject *data_generator,
                           bool write_legacy_png)
{
  png_structp png_ptr = NULL;
  png_infop info_ptr = NULL;
  PyObject * result = NULL;
  int bpc;
  FILE * fp = NULL;
  PyObject *iterator = NULL;

  /* TODO: try if this silliness helps
#if defined(PNG_LIBPNG_VER) && (PNG_LIBPNG_VER >= 10200)
  png_uint_32 mask, flags;
  
  flags = png_get_asm_flags(png_ptr);
  mask = png_get_asm_flagmask(PNG_SELECT_READ | PNG_SELECT_WRITE);
  png_set_asm_flags(png_ptr, flags | mask);
#endif
  */

  bpc = 8;
  
  fp = fopen(filename, "wb");
  if (!fp) {
    PyErr_SetFromErrno(PyExc_IOError);
    //PyErr_Format(PyExc_IOError, "Could not open PNG file for writing: %s", filename);
    goto cleanup;
  }

  png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, (png_voidp)NULL, png_write_error_callback, NULL);
  if (!png_ptr) {
    PyErr_SetString(PyExc_MemoryError, "png_create_write_struct() failed");
    goto cleanup;
  }

  info_ptr = png_create_info_struct(png_ptr);
  if (!info_ptr) {
    PyErr_SetString(PyExc_MemoryError, "png_create_info_struct() failed");
    goto cleanup;
  }

  if (setjmp(png_jmpbuf(png_ptr))) {
    goto cleanup;
  }

  png_init_io(png_ptr, fp);

  png_set_IHDR (png_ptr, info_ptr,
                w, h, bpc,
                has_alpha ? PNG_COLOR_TYPE_RGB_ALPHA : PNG_COLOR_TYPE_RGB,
                PNG_INTERLACE_NONE,
                PNG_COMPRESSION_TYPE_BASE,
                PNG_FILTER_TYPE_BASE);

  if (! write_legacy_png) {
    // Internal data is sRGB by the time it gets here.
    // Explicitly save with the recommended chunks to advertise that fact.
    png_set_sRGB_gAMA_and_cHRM (png_ptr, info_ptr, PNG_sRGB_INTENT_PERCEPTUAL);
  }

  // default (all filters enabled):                 1350ms, 3.4MB
  //png_set_filter(png_ptr, 0, PNG_FILTER_NONE);  // 790ms, 3.8MB
  //png_set_filter(png_ptr, 0, PNG_FILTER_PAETH); // 980ms, 3.5MB
  png_set_filter(png_ptr, 0, PNG_FILTER_SUB);     // 760ms, 3.4MB

  //png_set_compression_level(png_ptr, 0); // 0.49s, 32MB
  //png_set_compression_level(png_ptr, 1); // 0.98s, 9.6MB
  png_set_compression_level(png_ptr, 2);   // 1.08s, 9.4MB
  //png_set_compression_level(png_ptr, 9); // 18.6s, 9.3MB

  png_write_info(png_ptr, info_ptr);

  if (!has_alpha) {
    // input array format format is rgbu
    png_set_filler(png_ptr, 0, PNG_FILLER_AFTER);
  }

  {
    iterator = PyObject_GetIter(data_generator);
    if (!iterator) goto cleanup;

    int y = 0;
    while (y < h) {
      int rows;
      PyObject * arr = PyIter_Next(iterator);
      if (PyErr_Occurred()) goto cleanup;
      assert(arr); // iterator should have data
      assert(PyArray_ISALIGNED(arr));
      assert(PyArray_NDIM(arr) == 3);
      assert(PyArray_DIM(arr, 1) == w);
      assert(PyArray_DIM(arr, 2) == 4); // rgbu
      assert(PyArray_TYPE(arr) == NPY_UINT8);
      assert(PyArray_STRIDE(arr, 1) == 4);
      assert(PyArray_STRIDE(arr, 2) == 1);

      rows = PyArray_DIM(arr, 0);
//.........这里部分代码省略.........

static PyObject *somap(PyObject *self, PyObject *args)
{
    int ninputs, nsteps, nrows, ncols, depth, width, height, i, j, k;
    PyArrayObject *inputs_obj, *spread_obj;
    PyObject *somap_obj;
    double **inputs, **spread, ***somap;
    time_t t0, t1;
    npy_intp dims[3];

    if (!PyArg_ParseTuple(args, "O!O!iii", &PyArray_Type, &inputs_obj, &PyArray_Type, &spread_obj, &nrows, &ncols, &nsteps))
	RETURN_ERROR("need five arguments: inputs, spread, nrows, ncols, nsteps");

    if (!PyArray_Check(inputs_obj))
	RETURN_ERROR("inputs needs to be NumPy array");

    if (PyArray_NDIM(inputs_obj) != 2)
	RETURN_ERROR("inputs needs to be NumPy array with ndim 2");

    if (!PyArray_Check(spread_obj))
	RETURN_ERROR("spread needs to be NumPy array");

    if (PyArray_NDIM(spread_obj) != 2)
	RETURN_ERROR("spread needs to be NumPy array with ndim 2");

    if (PyArray_TYPE(inputs_obj) != NPY_DOUBLE)
	RETURN_ERROR("inputs needs to be array of doubles");

    if (PyArray_TYPE(spread_obj) != NPY_DOUBLE)
	RETURN_ERROR("spread needs to be array of doubles");

    ninputs = PyArray_DIM(inputs_obj, 0);
    depth = PyArray_DIM(inputs_obj, 1);

    width = PyArray_DIM(spread_obj, 0);
    height = PyArray_DIM(spread_obj, 1);

    //if (PyArray_AsCArray((PyObject**) &inputs_obj, (void *) &inputs, PyArray_DIMS(inputs_obj), PyArray_NDIM(inputs_obj), PyArray_DescrFromType(NPY_DOUBLE)) < 0)
    if (!(inputs = copyArray2(inputs_obj)))
	RETURN_ERROR("getting inputs as C array");

    //if (PyArray_AsCArray((PyObject**) &spread_obj, (void *) &spread, PyArray_DIMS(spread_obj), PyArray_NDIM(spread_obj), PyArray_DescrFromType(NPY_DOUBLE)) < 0)
    if (!(spread = copyArray2(spread_obj)))
    {
    	//PyArray_Free((PyObject*) inputs_obj, (void *) inputs);
        freeArray2(inputs, ninputs);
	RETURN_ERROR("getting spread as C array");
    }

    t0 = time(NULL);

    somap = selfOrganisingMap(inputs, ninputs, depth, nrows, ncols, spread, width, height, nsteps);

    t1 = time(NULL);

    printf("Time for just somap = %ld\n", t1-t0);

    //somap_obj = PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(NPY_DOUBLE), PyArray_NDIM(inputs_obj), PyArray_DIMS(inputs_obj), NULL, (void *) somap, 0, NULL);
    // below does not work because data not contiguous
    //somap_obj = PyArray_SimpleNewFromData(PyArray_NDIM(inputs_obj), PyArray_DIMS(inputs_obj), NPY_DOUBLE, (void *) somap);
    dims[0] = nrows;
    dims[1] = ncols;
    dims[2] = depth;
    somap_obj = PyArray_SimpleNew(3, dims, NPY_DOUBLE);
    for (i = 0; i < nrows; i++)
      for (j = 0; j < ncols; j++)
        for (k = 0; k < depth; k++)
          *((double *) PyArray_GETPTR3(somap_obj, i, j, k)) = somap[i][j][k];

    //PyArray_Free((PyObject*) inputs_obj, (void *) inputs);
    //PyArray_Free((PyObject*) spread_obj, (void *) spread);
    freeArray3(somap, nrows, ncols);
    freeArray2(inputs, ninputs);
    freeArray2(spread, width);

    return somap_obj;
}