From 23abdd884f76db602217b86b00b7ee0f3c045f8b Mon Sep 17 00:00:00 2001 From: Andrei Semena Date: Tue, 22 Jul 2025 13:42:28 +0300 Subject: [PATCH] found error in the estimation of the energy index for energy interpolation; added script to estimate likelihood for specified location --- chan_psf.c | 161 ++++++++++++++++++++++++++++++++++++++----- lkl_solver.py | 101 +++++++++++++++++++++++++++ source_detection2.py | 14 ++-- 3 files changed, 252 insertions(+), 24 deletions(-) create mode 100644 lkl_solver.py diff --git a/chan_psf.c b/chan_psf.c index 27b7590..c34cd19 100644 --- a/chan_psf.c +++ b/chan_psf.c @@ -180,21 +180,6 @@ static PyObject * solve_for_locations(PyObject *self, PyObject *args) 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; @@ -251,8 +236,6 @@ static PyObject * solve_for_locations(PyObject *self, PyObject *args) 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; }; @@ -347,7 +330,7 @@ static PyObject * solve_for_locations_eintp(PyObject *self, PyObject *args) long ctr, msum=0; double lkl, erf; - double pval, eloc, p2, p3, ptmp; + double pval, eloc, p2, p3; int idx1d, idx2d; //return psfdata + ((k*dims[1] + ei)*dims[2] + xi)*dims[3] + yi; @@ -462,11 +445,153 @@ static PyObject * solve_for_locations_eintp(PyObject *self, PyObject *args) +static PyObject * solve_for_rates(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, *eidx, *x, *y, *rates, *roll, *pk, *smat; + double xc, yc; + int rid, ctr; + double x1, y1, dx, dy, eloc; + + + if (!PyArg_ParseTuple(args, "OOOOOOOdddO", &psfi, &eidx, &x, &y, &roll, &pk, &rates, &xc, &yc, &eloc, &smat)) return NULL; + // -------------------------- =============== + // those are events properties those for sky smat it array for psf matrices + + npy_intp snew = {rates->dimensions[0]}; + PyArrayObject * lkls = PyArray_SimpleNew(1, &snew, NPY_DOUBLE); + double * lklsd = (double*) lkls->data; + + double * smatd = (double*) smat->data; + + 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, rate; + double * pkd = (double*) pk->data; + long * k = (long*)psfi->data; + double * ek = (double*)eidx->data; + int ei; + double* xptr = (double*) x->data; + double* yptr = (double*) y->data; + + + long msum=0; + double lkl, erf; + + double pval, p2, p3; + int idx1d, idx2d; + + for (rid=0; rid < rates->dimensions[0]; rid++) // loop over sky locations + { + msum = 0; + for (ctr=0; ctr < psfi->dimensions[0]; ctr++) // for each sky location loop over all provided events + { + x1 = (xc - xptr[ctr]); + y1 = (yc - 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 + ei = (int)(ek[ctr]); + erf = ek[ctr] - (double)(ei); + //printf("evt %d ei %d erf %f ek %f dx %f dy %f\n", ctr, ei, erf, ek[ctr], dx, dy); + + //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)); + + pval = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d, idx2d))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d, idx2d))*erf; + + //naive interpolation block + //------------------------------------------------------------------------------------------------------- + inpixdx = dx - (idx1d - 50); + inpixdy = dy - (idx2d - 50); + if (inpixdx > 0.) + { + p2 = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d + 1, idx2d))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d + 1, idx2d))*erf; + if (inpixdy > 0.) + { + p3 = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d, idx2d + 1))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d, idx2d + 1))*erf; + }else{ + inpixdy = -inpixdy; + p3 = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d, idx2d - 1))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d, idx2d - 1))*erf; + } + }else{ + p2 = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d - 1, idx2d))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d - 1, idx2d))*erf; + inpixdx = -inpixdx; + if (inpixdy > 0.) + { + p3 = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d, idx2d + 1))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d, idx2d + 1))*erf; + }else{ + inpixdy = -inpixdy; + p3 = (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei, idx1d, idx2d - 1))*(1. - erf) + (* eepsfvalfromptr(smatd, smat->dimensions, *(k + ctr), ei + 1, idx1d, idx2d - 1))*erf; + } + } + //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 + //------------------------------------------------------------------------------------------------------- + + if (pval > 1e-10) + { + bw[msum] = pval; + msum += 1; + //printf("%f %d\n", pval, msum); + }; + + }; + }; + + }; + if (msum > 0) + { + rate = (double) *((double*) rates->data + rid); + lkl = 0.; + for (ctr=0; ctr < msum; ctr ++) + { + lkl = lkl + log(rate*bw[ctr] + 1.); + } + *(lklsd + rid) = lkl - eloc*rate; + }else{ + *(lklsd + rid) = 0.; + }; + }; + //printf("loop done\n"); + + Py_END_ALLOW_THREADS; + + free(bw); + + PyObject *res = Py_BuildValue("O", lkls); + Py_DECREF(lkls); + return res; +} + + static PyMethodDef PSFMethods[] = { {"solve_for_locations", solve_for_locations, METH_VARARGS, "get coordinates within pixel based on its coordinates"}, {"solve_for_locations_eintp", solve_for_locations_eintp, METH_VARARGS, "compute likelihood using psf energy interpolation"}, {"put_psf_on_img", put_psf_on, METH_VARARGS, "put psf as is on img for all cooreindates "}, + {"solve_for_rates", solve_for_rates, METH_VARARGS, "computed likelihood at specified position for a series of rates"}, {NULL, NULL, 0, NULL} }; diff --git a/lkl_solver.py b/lkl_solver.py new file mode 100644 index 0000000..8f990e2 --- /dev/null +++ b/lkl_solver.py @@ -0,0 +1,101 @@ +import numpy as np +from astropy.io import fits +import matplotlib.pyplot as plt +import pickle +from astropy.wcs import WCS +import tqdm +from multiprocessing.pool import ThreadPool +from chan_psf import solve_for_locations, solve_for_locations_eintp, solve_for_rates + +psfe = np.array([1.8, 1.9, 3.0, 4.0, 6.0, 7.0, 8.0, 9.0]) + +def prepare_psf(evt): + """ + find all unique psf for observation and load in single 3d data cuve + return data cube with events slices indexes + """ + u, ui = np.unique(evt["psf_cube"], return_inverse=True) + data = np.array([np.load(p[3:])[:, ::-1,::-1].copy() for p in u]) + return data, ui + +def select_xychunksize(wcs, halfpsfsize=36./3600.): + """ + get wcs and find wcs pixel size of psf reach + """ + sizex = int(abs(halfpsfsize/wcs.wcs.cdelt[1])) + 2 + sizey = int(abs(halfpsfsize/wcs.wcs.cdelt[0])) + 2 + return sizex, sizey + +def read_wcs(h): + """ + read events wcs header + """ + w = WCS(naxis=2) + w.wcs.ctype = [h["TCTYP11"], h["TCTYP12"]] + w.wcs.crval = [h["TCRVL11"], h["TCRVL12"]] + w.wcs.cdelt = [h["TCDLT11"], h["TCDLT12"]] + w.wcs.crpix = [h["TCRPX11"], h["TCRPX12"]] + w = WCS(w.to_header()) + return w + + +def create_neighboring_blocks(x, y, sizex, sizey): + """ + schematically all sky is splitted on squares, which are approximatelly ~ 10 times greater then the psf + events for corresponding square are joined :: squer + diluttaion of psf reach + + coordinate system with events and all required coefficiets are fed to psf solver + current psf size is 25*0.5 arcsec (with up to \sqrt(2) factor in case of worst rolls -> 36'' + + """ + """ + event list already contains x and y for each event + """ + iix = (x//sizex + 0.5).astype(int) + iiy = (y//sizey + 0.5).astype(int) + isx, isy = np.mgrid[-1:2:1, -1:2:1] + oidx = np.repeat(np.arange(x.size), 9) + xyu, iixy, xyc = np.unique(np.array([np.repeat(iix, 9) + np.tile(isx.ravel(), x.size), + np.repeat(iiy, 9)+ np.tile(isy.ravel(), x.size)]), axis=1, return_counts=True, return_inverse=True) + + sord = np.argsort(iixy) + return oidx[sord], xyu, xyc + + +def lkls_for_rates(evt, expv, wcs, srcx, srcy, rates): + sizex, sizey = select_xychunksize(wcs) + x, y = evt["x"].astype(float), evt["y"].astype(float) + mask = np.logical_and.reduce([x > srcx - sizex//2, y > srcy - sizey//2, x < srcx + sizex//2, y < srcy + sizey//2], axis=0) + print("mask sum", srcx, srcy, mask.sum()) + eloc = evt[mask] + pickle.dump(eloc, open("eloc.pkl", "wb")) + psfdata, ui = prepare_psf(eloc) + xe, ye = np.copy(x[mask]), np.copy(y[mask]) + eidx = np.maximum(np.searchsorted(psfe*1e3, eloc["ENERGY"]) - 1, 0) + ee = np.maximum((eloc["ENERGY"]/1000. - psfe[eidx])/(psfe[eidx + 1] - psfe[eidx]), 0.).astype(float) + eidx + pk = np.copy(eloc["src_spec"]/eloc["bkg_spec"]).astype(float) + roll = np.copy(np.deg2rad(eloc["roll_pnt"])).astype(float) + #"OOOOOOOdddO", &psfi, &eidx, &x, &y, &roll, &pk, &rates, &xc, &yc, &eloc, &smat + # O O O O O O O d d d O" + print(ui, ee, xe, ye, roll, pk) + lkls = solve_for_rates(ui, ee, xe, ye, roll, pk, rates, srcx, srcy, expv, psfdata) + return lkls + + +if __name__ == "__main__": + p1 = fits.open("test.fits") + ewcs = read_wcs(p1[1].header) + wcs = WCS(fits.getheader("eR_spec_asp_0.fits.gz", 0)) + xc, yc = 4290, 4147 + xc, yc = 4643, 4223.7 + #xc, yc = 4147,4290 + xc, yc = ewcs.all_world2pix(wcs.all_pix2world([[xc, yc],], 0), 0).T + print(xc, yc) + eloc = 0.025 #0.0283 + #rates = np.array([4.2/eloc,]) #np.logspace(-0.5, 0.5, 129)*4.2/eloc + rates = np.logspace(-0.5, 0.5, 129)*1352/eloc #*4.2/eloc + lkls = lkls_for_rates(p1[1].data, eloc, ewcs, xc, yc, rates) + plt.plot(rates, lkls) + plt.axvline(rates[64]) + plt.show() + diff --git a/source_detection2.py b/source_detection2.py index 64d8639..b31cd8c 100644 --- a/source_detection2.py +++ b/source_detection2.py @@ -71,10 +71,10 @@ def make_srccount_and_detmap(emap, evt, h, wcs=None): x, y = evt["x"], evt["y"] else: ewcs = read_wcs(h) - x, y = wcs.all_world2pix(ewcs.all_pix2world(np.array([x, y]).T, 0), 0).T + x, y = wcs.all_world2pix(ewcs.all_pix2world(np.array([evt["x"], evt["y"]]).T, 0), 0).T - eidx = np.searchsorted(psfe*1e3, evt["ENERGY"]) - eidx = np.maximum((evt["ENERGY"]/1000. - psfe[eidx])/(psfe[eidx + 1] - psfe[eidx]), 0.) + eidx = np.maximum(np.searchsorted(psfe*1e3, evt["ENERGY"]) - 1, 0) + eidx = np.maximum((evt["ENERGY"]/1000. - psfe[eidx])/(psfe[eidx + 1] - psfe[eidx]), 0.) + eidx sizex, sizey = select_xychunksize(wcs) iidx, xyu, cts = create_neighboring_blocks(x, y, sizex, sizey) cc = np.zeros(cts.size + 1, int) @@ -115,8 +115,10 @@ def make_srccount_and_detmap(emap, evt, h, wcs=None): if __name__ == "__main__": p1 = fits.open("test.fits") #emap = fits.getdata("exp.map.gz") #np.full((8192, 8192), 10000.) - emap = fits.getdata("eR_spec_asp_0.fits.gz") #np.full((8192, 8192), 10000.) + emapf = fits.open("eR_spec_asp_0.fits.gz") #np.full((8192, 8192), 10000.) + emap = emapf[0].data + w = WCS(emapf[0].header) - wcs, cmap, pmap = make_srccount_and_detmap(emap, p1[1].data, p1[1].header) - fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(pmap - cmap, header=p1[1].header), fits.ImageHDU(cmap, header=p1[1].header)]).writeto("tmap4.fits.gz", overwrite=True) + wcs, cmap, pmap = make_srccount_and_detmap(emap, p1[1].data, p1[1].header, wcs=w) + fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(pmap - cmap, header=p1[1].header), fits.ImageHDU(cmap, header=p1[1].header)]).writeto("tmap5.fits.gz", overwrite=True) #fits.ImageHDU(data=pmap, header=wcs.to_header()).writeto("tmap4.fits.gz", overwrite=True)