209 lines
8.8 KiB
C
209 lines
8.8 KiB
C
#include "numpy/arrayobject.h"
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <stdbool.h>
|
|
#include <rate_solvers.h>
|
|
#ifndef MATH
|
|
#define MATh
|
|
#include <math.h>
|
|
#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);
|
|
};
|