initial commit
This commit is contained in:
adsam.job
parent
ef7e0d3f11
commit
b202022277
3
READMY.md
Normal file
3
READMY.md
Normal file
@ -0,0 +1,3 @@
|
||||
likelihood ratio source detection is working based on the psf and background model of the telescopes.
|
||||
It compares likelihood of the data for two models -- empty sky and sky with a source.
|
||||
The empty sky should produce well understandable distribution of measurements, which allows to obtain roughy constant rate of false sources with varying sensitivity, thus maximizing the ratio of detected and false sources.
|
||||
208
chan_psf.c
Normal file
208
chan_psf.c
Normal file
@ -0,0 +1,208 @@
|
||||
#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);
|
||||
};
|
||||
300
rate_solvers.c
Normal file
300
rate_solvers.c
Normal file
@ -0,0 +1,300 @@
|
||||
|
||||
double lkl_rate_condition(double r, double * s, double * b, int size)
|
||||
{
|
||||
int i;
|
||||
double res;
|
||||
res = 0.;
|
||||
for (i = 0; i < size; i++)
|
||||
{
|
||||
res += s[i]/(s[i]*r + b[i]);
|
||||
};
|
||||
return res;
|
||||
}
|
||||
|
||||
double lkl_rate_condition_pkr(double r, double * pk, int size)
|
||||
{
|
||||
int i;
|
||||
double res;
|
||||
res = 0.;
|
||||
for (i = 0; i < size; i++)
|
||||
{
|
||||
res += pk[i]/(pk[i]*r + 1.);
|
||||
};
|
||||
return res;
|
||||
}
|
||||
|
||||
double get_lkl_pkr(double r, double *pk, int size)
|
||||
{
|
||||
int i;
|
||||
double res;
|
||||
res = 0.;
|
||||
for (i = 0; i < size; i++)
|
||||
{
|
||||
res = res + log(pk[i]*r + 1.);
|
||||
};
|
||||
return res;
|
||||
}
|
||||
|
||||
|
||||
|
||||
double r8_epsilon ( )
|
||||
|
||||
/****************************************************************************
|
||||
|
||||
Purpose:
|
||||
|
||||
R8_EPSILON returns the R8 roundoff unit.
|
||||
|
||||
Discussion:
|
||||
|
||||
The roundoff unit is a number R which is a power of 2 with the
|
||||
property that, to the precision of the computer's arithmetic,
|
||||
1 < 1 + R
|
||||
but
|
||||
1 = ( 1 + R / 2 )
|
||||
|
||||
Licensing:
|
||||
|
||||
This code is distributed under the GNU LGPL license.
|
||||
|
||||
Modified:
|
||||
|
||||
01 September 2012
|
||||
|
||||
Author:
|
||||
|
||||
John Burkardt
|
||||
|
||||
Parameters:
|
||||
|
||||
Output, double R8_EPSILON, the R8 round-off unit.
|
||||
*/
|
||||
{
|
||||
const double value = 2.220446049250313E-016;
|
||||
return value;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
double brent (double rate, double exposure, double *sp, double *bp, int ssize)
|
||||
//****************************************************************************
|
||||
//
|
||||
// Purpose:
|
||||
//
|
||||
// ZERO seeks the root of a function F(X) in an interval [A,B].
|
||||
//
|
||||
// Discussion:
|
||||
//
|
||||
// The interval [A,B] must be a change of sign interval for F.
|
||||
// That is, F(A) and F(B) must be of opposite signs. Then
|
||||
// assuming that F is continuous implies the existence of at least
|
||||
// one value C between A and B for which F(C) = 0.
|
||||
//
|
||||
// The location of the zero is determined to within an accuracy
|
||||
// of 6 * MACHEPS * abs ( C ) + 2 * T.
|
||||
//
|
||||
// Thanks to Thomas Secretin for pointing out a transcription error in the
|
||||
// setting of the value of P, 11 February 2013.
|
||||
//
|
||||
// Licensing:
|
||||
//
|
||||
// This code is distributed under the GNU LGPL license.
|
||||
//
|
||||
// Modified:
|
||||
//
|
||||
// 11 February 2013
|
||||
//
|
||||
// Author:
|
||||
//
|
||||
// Original FORTRAN77 version by Richard Brent.
|
||||
// C++ version by John Burkardt.
|
||||
//
|
||||
// Reference:
|
||||
//
|
||||
// Richard Brent,
|
||||
// Algorithms for Minimization Without Derivatives,
|
||||
// Dover, 2002,
|
||||
// ISBN: 0-486-41998-3,
|
||||
// LC: QA402.5.B74.
|
||||
//
|
||||
// Parameters:
|
||||
//
|
||||
// Input, double A, B, the endpoints of the change of sign interval.
|
||||
//
|
||||
// Input, double T, a positive error tolerance.
|
||||
//
|
||||
// Input, func_base& F, the name of a user-supplied c++ functor
|
||||
// whose zero is being sought. The input and output
|
||||
// of F() are of type double.
|
||||
//
|
||||
// Output, double ZERO, the estimated value of a zero of
|
||||
// the function F.
|
||||
//
|
||||
{
|
||||
double a, b, t;
|
||||
double c;
|
||||
double d;
|
||||
double e;
|
||||
double fa;
|
||||
double fb;
|
||||
double fc;
|
||||
double m;
|
||||
double macheps;
|
||||
double p;
|
||||
double q;
|
||||
double r;
|
||||
double s;
|
||||
double sa;
|
||||
double sb;
|
||||
double tol;
|
||||
int i = 0;
|
||||
//
|
||||
// Make local copies of A and B.
|
||||
//
|
||||
a = 0.;
|
||||
b = ssize/exposure;
|
||||
t = 1e-7;
|
||||
|
||||
sa = a;
|
||||
sb = b;
|
||||
fa = lkl_rate_condition(sa, sp, bp, ssize) - exposure;
|
||||
fb = lkl_rate_condition(sb, sp, bp, ssize) - exposure;
|
||||
|
||||
c = sa;
|
||||
fc = fa;
|
||||
e = sb - sa;
|
||||
d = e;
|
||||
|
||||
macheps = r8_epsilon ( );
|
||||
|
||||
for ( ; ; )
|
||||
{
|
||||
i += 1;
|
||||
if ( fabs ( fc ) < fabs ( fb ) )
|
||||
{
|
||||
sa = sb;
|
||||
sb = c;
|
||||
c = sa;
|
||||
fa = fb;
|
||||
fb = fc;
|
||||
fc = fa;
|
||||
}
|
||||
|
||||
tol = 2.0 * macheps * fabs ( sb ) + t;
|
||||
m = 0.5 * ( c - sb );
|
||||
|
||||
if ( fabs ( m ) <= tol || fb == 0.0 )
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
if ( fabs ( e ) < tol || fabs ( fa ) <= fabs ( fb ) )
|
||||
{
|
||||
e = m;
|
||||
d = e;
|
||||
}
|
||||
else
|
||||
{
|
||||
s = fb / fa;
|
||||
|
||||
if ( sa == c )
|
||||
{
|
||||
p = 2.0 * m * s;
|
||||
q = 1.0 - s;
|
||||
}
|
||||
else
|
||||
{
|
||||
q = fa / fc;
|
||||
r = fb / fc;
|
||||
p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
|
||||
q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
|
||||
}
|
||||
|
||||
if ( 0.0 < p )
|
||||
{
|
||||
q = - q;
|
||||
}
|
||||
else
|
||||
{
|
||||
p = - p;
|
||||
}
|
||||
|
||||
s = e;
|
||||
e = d;
|
||||
|
||||
if ( 2.0 * p < 3.0 * m * q - fabs ( tol * q ) &&
|
||||
p < fabs ( 0.5 * s * q ) )
|
||||
{
|
||||
d = p / q;
|
||||
}
|
||||
else
|
||||
{
|
||||
e = m;
|
||||
d = e;
|
||||
}
|
||||
}
|
||||
sa = sb;
|
||||
fa = fb;
|
||||
|
||||
if ( tol < fabs ( d ) )
|
||||
{
|
||||
sb = sb + d;
|
||||
}
|
||||
else if ( 0.0 < m )
|
||||
{
|
||||
sb = sb + tol;
|
||||
}
|
||||
else
|
||||
{
|
||||
sb = sb - tol;
|
||||
}
|
||||
|
||||
fb = lkl_rate_condition(sb, sp, bp, ssize) - exposure;
|
||||
|
||||
if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
|
||||
{
|
||||
c = sa;
|
||||
fc = fa;
|
||||
e = sb - sa;
|
||||
d = e;
|
||||
}
|
||||
}
|
||||
/*printf("it counter %d\n", i);*/
|
||||
return sb;
|
||||
}
|
||||
|
||||
double get_phc_solution(double r, double e, double *s, double *b, int size)
|
||||
{
|
||||
int i;
|
||||
double rnew = r;
|
||||
for (i = 0; i < 200; i++)
|
||||
{
|
||||
rnew = lkl_rate_condition(r, s, b, size)*r/e;
|
||||
if ((fabs(rnew - r) < 1e-7) | (rnew*e < 1.)) break;
|
||||
r = rnew;
|
||||
};
|
||||
/*printf("it counter %d\n", i);*/
|
||||
|
||||
//printf("%d %e %e\n", size, s[0], r);
|
||||
return rnew;
|
||||
}
|
||||
|
||||
double get_phc_solution_pkr(double r, double e, double *pk, int size)
|
||||
{
|
||||
int i;
|
||||
double rnew = r;
|
||||
for (i = 0; i < 200; i++)
|
||||
{
|
||||
rnew = lkl_rate_condition_pkr(r, pk, size)*r/e;
|
||||
//printf("rnew %d %f\n", i, rnew*e);
|
||||
if ((fabs(rnew - r) < 1e-7) | (rnew*e < 0.001)) break;
|
||||
//if ((fabs(rnew - r) < 1e-7) | (rnew*e < 1.)) break;
|
||||
r = rnew;
|
||||
};
|
||||
/*printf("it counter %d\n", i);*/
|
||||
|
||||
//printf("%d %e %e\n", size, s[0], r);
|
||||
return rnew;
|
||||
}
|
||||
|
||||
14
rate_solvers.h
Normal file
14
rate_solvers.h
Normal file
@ -0,0 +1,14 @@
|
||||
#ifndef MATH
|
||||
#define MATh
|
||||
#include <math.h>
|
||||
#endif
|
||||
|
||||
double brent (double rate, double exposure, double *sp, double *bp, int ssize);
|
||||
|
||||
double get_phc_solution(double r, double e, double *s, double *b, int size);
|
||||
|
||||
double get_phc_solution_pkr(double r, double e, double *pk, int size);
|
||||
|
||||
double get_lkl_pkr(double r, double *pk, int size);
|
||||
|
||||
typedef double (*fptr)(double, double, double *, double *, int);
|
||||
15
setup.py
Normal file
15
setup.py
Normal file
@ -0,0 +1,15 @@
|
||||
from distutils.core import setup, Extension
|
||||
import numpy
|
||||
|
||||
m1 = Extension("src_rate_solvers", ["rate_solvers.c", "src_rate_solvers.c"], include_dirs = ["./", numpy.get_include()])
|
||||
m2 = Extension("chan_psf", ["rate_solvers.c", "chan_psf.c"], include_dirs = ["./", numpy.get_include()])
|
||||
|
||||
def main():
|
||||
setup(name="lkl_solution",
|
||||
version="0.1",
|
||||
author="andrey",
|
||||
author_email="san@iki.rssi.ru",
|
||||
ext_modules=[m1, m2])
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
207
source_detection.py
Normal file
207
source_detection.py
Normal file
@ -0,0 +1,207 @@
|
||||
import asyncio
|
||||
import numpy as np
|
||||
from astropy.io import fits
|
||||
import matplotlib.pyplot as plt
|
||||
from astropy.wcs import WCS
|
||||
import tqdm
|
||||
from multiprocessing.pool import ThreadPool
|
||||
from chan_psf import solve_for_locations
|
||||
|
||||
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_path"], return_inverse=True)
|
||||
data = np.array([np.load(p)[::-1,::-1] 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
|
||||
print(sizex, sizey)
|
||||
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 make_srccount_and_detmap(emap, evt, h, wcs=None):
|
||||
psfdata, ui = prepare_psf(evt)
|
||||
if wcs is None:
|
||||
wcs = read_wcs(h)
|
||||
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
|
||||
|
||||
sizex, sizey = select_xychunksize(wcs)
|
||||
iidx, xyu, cts = create_neighboring_blocks(x, y, sizex, sizey)
|
||||
cc = np.zeros(cts.size + 1, int)
|
||||
cc[1:] = np.cumsum(cts)
|
||||
cmap, pmap = np.zeros(emap.shape, float), np.zeros(emap.shape, float)
|
||||
#xe, ye, pk, roll, psfi = np.copy(evt["x"][iidx]), np.copy(evt["y"][iidx]), np.copy((evt["quant_eff"]/evt["bkg_model"])[iidx]), np.copy(evt["roll_pnt"][iidx]), np.copy(ui[iidx])
|
||||
xe = np.copy(x[iidx]).astype(float)
|
||||
ye = np.copy(y[iidx]).astype(float)
|
||||
pk = np.copy((2e-3/evt["bkg_model"])[iidx]).astype(float)
|
||||
roll = np.copy(np.deg2rad(evt["roll_pnt"][iidx])).astype(float)
|
||||
psfi = np.copy(ui[iidx])
|
||||
|
||||
yg, xg = np.mgrid[0:sizey:1, 0:sizex:1]
|
||||
|
||||
def worker(ixys):
|
||||
i, (xs, ys) = ixys
|
||||
eloc = emap[ys*sizey:ys*sizey+sizey, xs*sizex:xs*sizex+sizex]
|
||||
mask = eloc > 0.
|
||||
xl = (xg[mask] + xs*sizex).astype(float)
|
||||
yl = (yg[mask] + ys*sizey).astype(float)
|
||||
ell = (eloc[mask]).astype(float)
|
||||
if np.any(mask):
|
||||
"""
|
||||
mask = np.any([(np.abs(x - xl) < 20) & (np.abs(y - yl) < 20) for x, y in zip(xe[cc[i]:cc[i+1]], ye[cc[i]:cc[i+1]])], axis=0)
|
||||
if mask.sum() == 0:
|
||||
return None
|
||||
xl, yl, ell = xl[mask], yl[mask], ell[mask]
|
||||
cr, pr = solve_for_locations(psfi[cc[i]:cc[i+1]], xe[cc[i]:cc[i+1]], ye[cc[i]:cc[i+1]], roll[cc[i]:cc[i+1]], pk[cc[i]:cc[i+1]], xl, yl, ell, psfdata)
|
||||
print(pr.min(), pr.max())
|
||||
plt.scatter(xl, yl, c=pr)
|
||||
plt.scatter(xe[cc[i]: cc[i + 1]], ye[cc[i]: cc[i+1]], marker="x", c="r", s=50, zorder=20)
|
||||
plt.show()
|
||||
"""
|
||||
cr, pr = solve_for_locations(psfi[cc[i]:cc[i+1]], xe[cc[i]:cc[i+1]], ye[cc[i]:cc[i+1]], roll[cc[i]:cc[i+1]], pk[cc[i]:cc[i+1]], xl, yl, ell, psfdata)
|
||||
else:
|
||||
xl, yl, cr, pr = np.empty(0, int),np.empty(0, int),np.empty(0, float),np.empty(0, float)
|
||||
return xl.astype(int), yl.astype(int), cr, pr
|
||||
|
||||
|
||||
"""
|
||||
for r in enumerate(xyu.T):
|
||||
worker(r)
|
||||
|
||||
"""
|
||||
tpool = ThreadPool(8)
|
||||
for xl, yl, cr, pr in tqdm.tqdm(tpool.imap_unordered(worker, enumerate(xyu.T)), total=xyu.shape[1]):
|
||||
cmap[yl, xl] = cr
|
||||
pmap[yl, xl] = pr
|
||||
|
||||
"""
|
||||
for i, (xs, ys) in enumerate(xyu.T):
|
||||
eloc = emap[xs:xs+sizex,ys:ys+sizey]
|
||||
mask = eloc > 0.
|
||||
xl = (xg[mask] + xs*sizex).astype(float)
|
||||
yl = (yg[mask] + ys*sizey).astype(float)
|
||||
ell = (eloc[mask]).astype(float)
|
||||
print(i, xs, ys, xl.min(), xl.max(), yl.min(), yl.max())
|
||||
cr, pr = solve_for_locations(psfi[cc[i]:cc[i+1]], xe[cc[i]:cc[i+1]], ye[cc[i]:cc[i+1]], roll[cc[i]:cc[i+1]], pk[cc[i]:cc[i+1]], xl, yl, ell, psfdata)
|
||||
xl, yl = xl.astype(int), yl.astype(int)
|
||||
cmap[xl, yl] = cr
|
||||
pmap[xl, yl] = pr
|
||||
i = 0
|
||||
xs, ys = xyu[:, i]
|
||||
eloc = emap[xs:xs+sizex,ys:ys+sizey]
|
||||
mask = eloc > 0.
|
||||
xl = (xg[mask] + xs*sizex).astype(float)
|
||||
yl = (yg[mask] + ys*sizey).astype(float)
|
||||
ell = (eloc[mask]).astype(float)
|
||||
print(i, xs, ys, xl.min(), xl.max(), yl.min(), yl.max())
|
||||
cr, pr = solve_for_locations(psfi[cc[i]:cc[i+1]], xe[cc[i]:cc[i+1]], ye[cc[i]:cc[i+1]], roll[cc[i]:cc[i+1]], pk[cc[i]:cc[i+1]], xl, yl, ell, psfdata)
|
||||
xl, yl = xl.astype(int), yl.astype(int)
|
||||
cmap[xl, yl] = cr
|
||||
pmap[xl, yl] = pr
|
||||
"""
|
||||
return wcs, cmap, pmap
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
p1 = fits.open("updated_test.fits")
|
||||
emap = fits.getdata("exp.fits.gz") #np.full((8192, 8192), 10000.)
|
||||
|
||||
wcs, cmap, pmap = make_srccount_and_detmap(emap, p1[1].data, p1[1].header)
|
||||
fits.ImageHDU(data=pmap, header=wcs.to_header()).writeto("tmap3.fits.gz", overwrite=True)
|
||||
"""
|
||||
|
||||
if __name__ == "__main__":
|
||||
p1 = fits.open("updated_test.fits")
|
||||
evt = p1[1].data
|
||||
|
||||
emap = np.full((8192, 8192), 14139.9)
|
||||
psfdata, ui = prepare_psf(evt)
|
||||
#psfdata = np.zeros(psfdata.shape, float)
|
||||
#psfdata[:, 50, 50] = 1.
|
||||
wcs = read_wcs(p1[1].header)
|
||||
|
||||
sizex, sizey = select_xychunksize(wcs)
|
||||
iidx, xyu, cts = create_neighboring_blocks(evt["x"], evt["y"], sizex, sizey)
|
||||
print(xyu.shape)
|
||||
cc = np.zeros(cts.size + 1, int)
|
||||
cc[1:] = np.cumsum(cts)
|
||||
cmap, pmap = np.zeros(emap.shape, float), np.zeros(emap.shape, float)
|
||||
#xe, ye, pk, roll, psfi = np.copy(evt["x"][iidx]), np.copy(evt["y"][iidx]), np.copy((evt["quant_eff"]/evt["bkg_model"])[iidx]), np.copy(evt["roll_pnt"][iidx]), np.copy(ui[iidx])
|
||||
xe, ye, pk, roll, psfi = np.copy(evt["x"][iidx]), np.copy(evt["y"][iidx]), np.copy((2e-3/evt["bkg_model"])[iidx]), np.copy(np.deg2rad(evt["roll_pnt"][iidx])), np.copy(ui[iidx])
|
||||
xe = xe.astype(float)
|
||||
ye = ye.astype(float)
|
||||
roll = roll.astype(float)
|
||||
pk = pk.astype(float)
|
||||
print("pk.dtype", pk.dtype, "xe.dtype", xe.dtype)
|
||||
|
||||
#plt.hist(pk, 100)
|
||||
#plt.show()
|
||||
xg, yg = np.mgrid[0:sizex:1, 0:sizey:1]
|
||||
|
||||
i = 60# np.argmax(cts)
|
||||
xs, ys = xyu[:, i]
|
||||
|
||||
eloc = emap[xs:xs+sizex,ys:ys+sizey]
|
||||
mask = eloc > 0.
|
||||
|
||||
xl = (xg[mask] + xs*sizex).astype(float)
|
||||
yl = (yg[mask] + ys*sizey).astype(float)
|
||||
plt.scatter(xe[cc[i]:cc[i+1]], ye[cc[i]: cc[i + 1]], marker="+", color="r")
|
||||
plt.scatter(xl, yl, marker="x", color="k", zorder=20)
|
||||
plt.show()
|
||||
|
||||
ell = np.copy(eloc[mask])
|
||||
print("check pk in python", pk[[cc[i], cc[i] + 10, cc[i] + 20]])
|
||||
print("check pk in python", psfi[[cc[i], cc[i] + 10, cc[i] + 20]])
|
||||
print("check pk in python", xe[[cc[i], cc[i] + 10, cc[i] + 20]])
|
||||
cr, pr = solve_for_locations(psfi[cc[i]:cc[i+1]], xe[cc[i]:cc[i+1]], ye[cc[i]:cc[i+1]], roll[cc[i]:cc[i+1]], pk[cc[i]:cc[i+1]], xl, yl, ell, psfdata)
|
||||
plt.scatter(xl, yl, c=pr)
|
||||
plt.show()
|
||||
|
||||
|
||||
"""
|
||||
8
src_rate_solver.h
Normal file
8
src_rate_solver.h
Normal file
@ -0,0 +1,8 @@
|
||||
|
||||
double brent (double rate, double exposure, double *sp, double *bp, int ssize);
|
||||
|
||||
double get_phc_solution(double r, double e, double *s, double *b, int size);
|
||||
|
||||
double get_phc_solution_pkr(double r, double e, double *pk, int size);
|
||||
|
||||
double get_lkl_pkr(double r, double *pk, int size);
|
||||
Loading…
x
Reference in New Issue
Block a user