found error in the estimation of the energy index for energy interpolation; added script to estimate likelihood for specified location

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}
 }; 
 

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()
+

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(np.searchsorted(psfe*1e3, evt["ENERGY"]) - 1, 0)
-    eidx = np.maximum((evt["ENERGY"]/1000. - psfe[eidx])/(psfe[eidx + 1] - psfe[eidx]), 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)
+    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("tmap4.fits.gz", overwrite=True)
+    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)