#include "numpy/arrayobject.h" #include #include #include #include #ifndef MATH #define MATh #include #endif double * psfvalfromptr(double * psfdata, npy_intp * dims, int k, int xi, int yi) { //printf("check didx %d %.2e\n", ((dims[1]*eidx + k)*dims[2] + xi)*dims[3] + yi, *( psfdata + ((dims[1]*eidx + k)*dims[2] + xi)*dims[3] + yi)); return psfdata + (k*dims[1] + xi)*dims[2] + yi; }; static PyObject * solve_for_locations(PyObject *self, PyObject *args) { //xc, yc --- wcs locations, events has coordinates in the same locations, and psf have the same grid as well // the only additional parameter to events are pk scale (rate scale in respect to psf) and rotation angle PyArrayObject *psfi, *x, *y, *roll, *pk, *xc, *yc, *smat, *emap; int loc, ctr; double x1, y1, dx, dy; if (!PyArg_ParseTuple(args, "OOOOOOOOO", &psfi, &x, &y, &roll, &pk, &xc, &yc, &emap, &smat)) return NULL; // -------------------------- =============== // those are events properties those for sky smat it array for psf matrices npy_intp snew = {xc->dimensions[0]}; PyArrayObject * cmap = PyArray_SimpleNew(1, &snew, NPY_DOUBLE); PyArrayObject * pmap = PyArray_SimpleNew(1, &snew, NPY_DOUBLE); double * cmapd = (double*) cmap->data; double * pmapd = (double*) pmap->data; double * smatd = (double*) smat->data; //printf("smat 1d %d xc size %d x size %d\n", smat->dimensions[0], xc->dimensions[0], x->dimensions[0]); double *ca = (double*)malloc(sizeof(double)*x->dimensions[0]); double *sa = (double*)malloc(sizeof(double)*x->dimensions[0]); double* nparrptr = (double*) roll->data; for (ctr=0; ctr < x->dimensions[0]; ctr++) { ca[ctr] = cos(nparrptr[ctr]); sa[ctr] = sin(nparrptr[ctr]); }; double * bw = (double*)malloc(sizeof(double)*psfi->dimensions[0]); //not more then thet will be used for each location Py_BEGIN_ALLOW_THREADS; double inpixdx, inpixdy; double * pkd = (double*) pk->data; long * k = (long*)psfi->data; double* xptr = (double*) x->data; double* yptr = (double*) y->data; double *xcptr = (double*) xc->data; double *ycptr = (double*) yc->data; long ctr, msum=0; double lkl; double pval, eloc, p2, p3; int idx1d, idx2d; /* pval = *psfvalfromptr(smatd, smat->dimensions, 1744, 50, 50); printf("1744 50 50 %f\n", pval); pval = *psfvalfromptr(smatd, smat->dimensions, 1744, 40, 48); printf("1744 40 48 %f\n", pval); pval = *psfvalfromptr(smatd, smat->dimensions, 1744, 20, 52); printf("1744 20 52 %f\n", pval); printf("bwd %f %f %f\n", xptr[0], xptr[10], xptr[20]); printf("xptr %f %f %f\n",xcptr[0], xcptr[10], xcptr[20]); printf("dimension %d\n", xc->dimensions[0]); */ for (loc=0; loc < xc->dimensions[0]; loc++) // loop over sky locations { msum = 0; for (ctr=0; ctr < psfi->dimensions[0]; ctr++) // for each sky location loop over all provided events { x1 = (xcptr[loc] - xptr[ctr]); y1 = (ycptr[loc] - yptr[ctr]); //rotate by the event roll angle, dx dy centered at the psf center (central pixel of 101x101 map) dx = x1*ca[ctr] - y1*sa[ctr]; //+ 50; dy = y1*ca[ctr] + x1*sa[ctr]; // + 50.; // temporary hardcode psf shape is 101x101 //current psf shape is 101: if ((dx > -50) && (dx < 50)) { if ((dy > -50) && (dy < 50)) { idx1d = (int)((dx + 50.5)); // float dx from -0.5 to 0.5 should fell in the 50-th pixel idx2d = (int)((dy + 50.5)); //printf("dx dy %f %f\n", dx, dy); pval = * psfvalfromptr(smatd, smat->dimensions, *(k + ctr), idx1d, idx2d); //naive interpolation block //------------------------------------------------------------------------------------------------------- inpixdx = dx - (idx1d - 50); inpixdy = dy - (idx2d - 50); if (inpixdx > 0.) { p2 = * psfvalfromptr(smatd, smat->dimensions, *(k + ctr), idx1d + 1, idx2d); if (inpixdy > 0.) { p3 = * psfvalfromptr(smatd, smat->dimensions, * (k + ctr), idx1d, idx2d + 1); }else{ inpixdy = -inpixdy; p3 = * psfvalfromptr(smatd, smat->dimensions, * (k + ctr), idx1d, idx2d - 1); } }else{ p2 = * psfvalfromptr(smatd, smat->dimensions, * (k + ctr), idx1d - 1, idx2d); inpixdx = -inpixdx; if (inpixdy > 0.) { p3 = * psfvalfromptr(smatd, smat->dimensions, * (k + ctr), idx1d, idx2d + 1); }else{ inpixdy = -inpixdy; p3 = * psfvalfromptr(smatd, smat->dimensions, * (k + ctr), idx1d, idx2d - 1); } } //printf("pval %f %f %f %f %f %d %d %d\n", pval, p2, p3, inpixdx, inpixdy, idx1d, idx2d, k[ctr]); pval = (pval + inpixdx*(p2 - pval) + inpixdy*(p3 - pval))* (*(pkd + ctr)); // interpolation up to here //------------------------------------------------------------------------------------------------------- //pval = pval * (*(pkd + ctr)); if (pval > 1e-10) { bw[msum] = pval; //printf("%d %d %d %f %f %f %f %f\n", k[ctr], idx1d, idx2d, inpixdx, inpixdy, dx, dy, pval); //printf("%d %d %d %f %f %f %f %f\n", k[ctr], idx1d, idx2d, xptr[ctr], yptr[ctr], xcptr[loc], ycptr[loc], pval); msum += 1; }; }; }; }; if (msum > 0) { eloc = (double) *((double*) emap->data + loc); pval = get_phc_solution_pkr((double) msum, eloc, bw, msum); *(cmapd + loc) = pval*eloc; lkl = 0.; for (ctr=0; ctr < msum; ctr ++) { lkl = lkl + log(pval*bw[ctr] + 1.); } //printf("loc %d %d %f %f %f %f\n", loc, msum, bw[0], eloc, pval, lkl); *(pmapd + loc) = lkl; //log(lkl); //get_lkl_pkr(pval, bw, msum); /* lkl = 0; for (ctr=0; ctr < msum; ctr ++) { lkl = lkl + bw[ctr]; } *(pmapd + loc) = lkl; printf("%d %d %f\n", msum, loc, pmapd[loc]); */ }else{ *(cmapd + loc) = 0.; *(pmapd + loc) = 0.; }; }; //printf("loop done\n"); Py_END_ALLOW_THREADS; free(bw); PyObject *res = Py_BuildValue("OO", cmap, pmap); Py_DECREF(cmap); Py_DECREF(pmap); return res; } static PyMethodDef PSFMethods[] = { {"solve_for_locations", solve_for_locations, METH_VARARGS, "get coordinates within pixel based on its coordinates"}, {NULL, NULL, 0, NULL} }; static struct PyModuleDef psf_c_module = { PyModuleDef_HEAD_INIT, "chan_psf", NULL, -1, PSFMethods }; PyMODINIT_FUNC PyInit_chan_psf(void) { assert(! PyErr_Occurred()); if (PyErr_Occurred()) {return NULL;} import_array(); return PyModule_Create(&psf_c_module); };