p61743098
/
sagecal
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

update

This commit is contained in:
Sarod Yatawatta 2017-08-16 12:32:05 +02:00
parent eb9769d462
commit 87d9f47f6e
10 changed files with 200 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/buildsky/annotate.py
View File

@ -44,7 +44,7 @@ def annotate_lsm_sky(infilename,clusterfilename,outfilename,clid=None,color='yel
\s+ # skip white space
(?P<col7>[-+]?\d+(\.\d+)?) # Dec angle - sec
\s+ # skip white space
(?P<col8>[-+]?\d+(\.\d+)?) # Stokes I - Flux
(?P<col8>[-+]?(\d+(\.\d*)?|\d*\.\d+)([eE][-+]?\d+)?) # Stokes I - Flux
\s+ # skip white space
(?P<col9>[-+]?(\d+(\.\d*)?|\d*\.\d+)([eE][-+]?\d+)?) # Stokes Q - Flux
\s+ # skip white space
@ -79,7 +79,7 @@ def annotate_lsm_sky(infilename,clusterfilename,outfilename,clid=None,color='yel
\s+ # skip white space
(?P<col7>[-+]?\d+(\.\d+)?) # Dec angle - sec
\s+ # skip white space
(?P<col8>[-+]?\d+(\.\d+)?) # Stokes I - Flux
(?P<col8>[-+]?(\d+(\.\d*)?|\d*\.\d+)([eE][-+]?\d+)?) # Stokes I - Flux
\s+ # skip white space
(?P<col9>[-+]?(\d+(\.\d*)?|\d*\.\d+)([eE][-+]?\d+)?) # Stokes Q - Flux
\s+ # skip white space
@ -111,17 +111,35 @@ def annotate_lsm_sky(infilename,clusterfilename,outfilename,clid=None,color='yel
SR={} # sources
for eachline in all:
v=pp.search(eachline)
v=pp1.search(eachline)
if v!= None:
# find RA,DEC (rad) and flux
mra=(float(v.group('col2'))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0
mdec=(float(v.group('col5'))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)
# find RA,DEC (rad) and flux, with proper sign
if (float(v.group('col2'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mra=mysign*(abs(float(v.group('col2')))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0
if (float(v.group('col5'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mdec=mysign*(abs(float(v.group('col5')))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)
SR[str(v.group('col1'))]=(mra,mdec,float(v.group('col8')))
else:
v=pp1.search(eachline)
v=pp.search(eachline)
if v!= None:
mra=(float(v.group('col2'))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0
mdec=(float(v.group('col5'))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)
if (float(v.group('col2'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mra=mysign*(abs(float(v.group('col2')))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0
if (float(v.group('col5'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mdec=mysign*(abs(float(v.group('col5')))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)
SR[str(v.group('col1'))]=(mra,mdec,float(v.group('col8')))
print 'Read %d sources'%len(SR)
@ -140,8 +158,8 @@ def annotate_lsm_sky(infilename,clusterfilename,outfilename,clid=None,color='yel
for eachline in all:
v=pp.search(eachline)
if v!= None:
# iterate over list of source names
CL[str(v.group('col1'))]=re.split('\W+',re.sub('\n','',str(v.group('col3'))))
# iterate over list of source names (names can also have a '.')
CL[str(v.group('col1'))]=re.split('[^a-zA-Z0-9_\.]+',re.sub('\n','',str(v.group('col3'))))
print 'Read %d clusters'%len(CL)

37
src/buildsky/buildmultisky.c
View File

@ -195,9 +195,10 @@ fillup_pixel_hashtablef(long totalrows, long offset, long firstrow, long nrows,
bmaj (rad),bmin (rad) ,pa (rad): PSF parameters from the FITS files: Nfx1 arrays
minpix: footprint of MEAN psf in pixels
ignlist: list of islands to ignore (integer list)
donegative: fit -ve pixels instead of positive
*/
int
read_fits_file_f(const char *fitsdir, const char *maskfile, GHashTable **pixtable, int *Nf, double **freqs, double **bmaj, double **bmin, double **pa, int beam_given, double *minpix, GList *ignlist) {
read_fits_file_f(const char *fitsdir, const char *maskfile, GHashTable **pixtable, int *Nf, double **freqs, double **bmaj, double **bmin, double **pa, int beam_given, double *minpix, GList *ignlist, int donegative) {
// hash table, list params
GHashTableIter iter;
GList *li;
@ -211,12 +212,10 @@ read_fits_file_f(const char *fitsdir, const char *maskfile, GHashTable **pixtabl
int status;
int naxis;
int bitpix;
long int new_naxis[4];
int ii,jj,kk;
int datatype=0;
long int totalpix;
double bscale,bzero;
long int increment[4]={1,1,1,1};
int null_flag=0;
@ -324,16 +323,6 @@ read_fits_file_f(const char *fitsdir, const char *maskfile, GHashTable **pixtabl
fbuff.arr_dims.hpix[1]=fbuff.arr_dims.d[1];
fbuff.arr_dims.hpix[2]=1; /* freq axes */
fbuff.arr_dims.hpix[3]=1;
/******* create new array **********/
new_naxis[0]=fbuff.arr_dims.hpix[0]-fbuff.arr_dims.lpix[0]+1;
new_naxis[1]=fbuff.arr_dims.hpix[1]-fbuff.arr_dims.lpix[1]+1;
new_naxis[2]=fbuff.arr_dims.hpix[2]-fbuff.arr_dims.lpix[2]+1;
new_naxis[3]=fbuff.arr_dims.hpix[3]-fbuff.arr_dims.lpix[3]+1;
/* calculate total number of pixels */
totalpix=((fbuff.arr_dims.hpix[0]-fbuff.arr_dims.lpix[0]+1)
*(fbuff.arr_dims.hpix[1]-fbuff.arr_dims.lpix[1]+1)
*(fbuff.arr_dims.hpix[2]-fbuff.arr_dims.lpix[2]+1)
*(fbuff.arr_dims.hpix[3]-fbuff.arr_dims.lpix[3]+1));
/* define input column structure members for the iterator function */
fits_iter_set_file(&cols[0], fbuff.fptr);
@ -654,10 +643,14 @@ printf("freq=%lf\n",(*freqs)[cnt]);
0, &mypix, &null_flag, &status);
/* if pixel has -ve flux, throw it away */
if (mypix<0.0) {
if (!donegative && mypix<0.0) {
mypix=0.0;
} else if (donegative && mypix >0.0) {
/* if pixel has +ve flux, throw it away */
mypix=0.0;
}
ppix->sI[cnt]=mypix;
/* always use +ve flux for fitting */
ppix->sI[cnt]=(donegative?-mypix:mypix);
val->stI+=mypix;
#ifdef DEBUG
printf("[pixel (%d,%d), lm (%lf,%lf), radec (%lf,%lf), sI ( ",ppix->x,ppix->y,ppix->l,ppix->m,ppix->ra,ppix->dec);
@ -904,6 +897,7 @@ add_guard_pixels_f(GList *pixlist, int Nf, double threshold, hpixelf **parr, int
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
return 1;
}
GList *pixval0=pixval;
ci=1;
#ifdef DEBUG
printf("Before x:");
@ -915,6 +909,7 @@ add_guard_pixels_f(GList *pixlist, int Nf, double threshold, hpixelf **parr, int
printf("%u ",*xkey);
#endif
}
g_list_free(pixval0);
#ifdef DEBUG
printf("\n");
#endif
@ -963,6 +958,7 @@ add_guard_pixels_f(GList *pixlist, int Nf, double threshold, hpixelf **parr, int
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
return 1;
}
pixval0=pixval;
ci=1;
#ifdef DEBUG
printf("Before y:");
@ -974,6 +970,7 @@ add_guard_pixels_f(GList *pixlist, int Nf, double threshold, hpixelf **parr, int
printf("%u ",*xkey);
#endif
}
g_list_free(pixval0);
#ifdef DEBUG
printf("\n");
#endif
@ -1276,7 +1273,7 @@ process_pixels_f(GHashTable *pixtable, int Nf, double *freqs, double *bmaj, doub
int
write_world_coords_f(const char *imgfile, GHashTable *pixtable, double minpix, int Nf, double *freqs, double *bmaj, double *bmin, double *bpa, double ref_freq, int outformat, double clusterratio, int nclusters, const char *unistr){
write_world_coords_f(const char *imgfile, GHashTable *pixtable, double minpix, int Nf, double *freqs, double *bmaj, double *bmin, double *bpa, double ref_freq, int outformat, double clusterratio, int nclusters, const char *unistr,int donegative, int scaleflux){
GHashTableIter iter;
GList *li,*pli;
pixellistf *val;
@ -1482,11 +1479,13 @@ write_world_coords_f(const char *imgfile, GHashTable *pixtable, double minpix, i
}
//printf("BEAM %lf,%lf\n",bmaj,bmin);
total_model_flux*=minpix; /* FIXME: no scale change */
fluxscale=val->stI/(total_model_flux);
if (scaleflux) {
fluxscale=val->stI/(total_model_flux);
}
}
printf("(%d) Total flux/beam=%lf model=%lf scale=%lf\n",*key_,val->stI/minpix,total_model_flux/minpix,fluxscale);
/* do not scale up flux */
fluxscale=1.0;
fluxscale=(donegative?-1.0:1.0);
/* try to cluster */
cluslist=NULL;
if (clusterratio>0.0) {
@ -1497,7 +1496,7 @@ write_world_coords_f(const char *imgfile, GHashTable *pixtable, double minpix, i
printf("Choosing clustered version ");
li=cluslist;
/* if cluster only have 1 source, normalize by peak flux, not the sum */
if (g_list_length(cluslist)==1) {
if (g_list_length(cluslist)==1 && scaleflux) {
srcx=li->data;
fluxscale=0.0;
for (ii=0;ii<Nf;++ii) {

12
src/buildsky/buildsky.c
View File

@ -559,7 +559,7 @@ read_fits_file(const char *imgfile, const char *maskfile, GHashTable **pixtable,
int
write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, double bmaj, double bmin, double bpa, int outformat, double clusterratio, int nclusters,const char *unistr){
write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, double bmaj, double bmin, double bpa, int outformat, double clusterratio, int nclusters,const char *unistr, int scaleflux){
GHashTableIter iter;
GList *li,*pli;
pixellist *val;
@ -722,7 +722,9 @@ write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, dou
}
//printf("BEAM %lf,%lf\n",bmaj,bmin);
total_model_flux*=1.1*minpix; /* FIXME: no scale change */
fluxscale=val->stI/(total_model_flux);
if (scaleflux) {
fluxscale=val->stI/(total_model_flux);
}
}
printf("(%d) Total flux/beam=%lf model=%lf scale=%lf\n",*key_,val->stI/minpix,total_model_flux/minpix,fluxscale);
/* try to cluster */
@ -736,7 +738,7 @@ write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, dou
printf("Choosing clustered version ");
li=cluslist;
/* if cluster only have 1 source, normalize by peak flux, not the sum */
if (g_list_length(cluslist)==1) {
if (g_list_length(cluslist)==1 && scaleflux) {
srcx=li->data;
/* peak absolute flux */
peak_abs=-1e6;
@ -1077,6 +1079,7 @@ add_guard_pixels(GList *pixlist, double threshold, hpixel **parr, int *n) {
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
return 1;
}
GList *pixval0=pixval;
ci=1;
#ifdef DEBUG
printf("Before x:");
@ -1088,6 +1091,7 @@ add_guard_pixels(GList *pixlist, double threshold, hpixel **parr, int *n) {
printf("%u ",*xkey);
#endif
}
g_list_free(pixval0);
#ifdef DEBUG
printf("\n");
#endif
@ -1136,6 +1140,7 @@ add_guard_pixels(GList *pixlist, double threshold, hpixel **parr, int *n) {
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
return 1;
}
pixval0=pixval;
ci=1;
#ifdef DEBUG
printf("Before y:");
@ -1147,6 +1152,7 @@ add_guard_pixels(GList *pixlist, double threshold, hpixel **parr, int *n) {
printf("%u ",*xkey);
#endif
}
g_list_free(pixval0);
#ifdef DEBUG
printf("\n");
#endif

10
src/buildsky/buildsky.h
View File

@ -285,9 +285,10 @@ filter_pixels(GHashTable *pixtable, double wcutoff);
clusterratio: components closer than clusteratio*(bmaj+bmin)/2 will be merged
nclusters: no of max clusters to cluster the sky
unistr: unique string for source names
scaleflux: if 1, scale model flux to match total flux of island
*/
extern int
write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, double bmaj, double bmin, double bpa, int outformat, double clusterratio, int nclusters,const char *unistr);
write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, double bmaj, double bmin, double bpa, int outformat, double clusterratio, int nclusters,const char *unistr, int scaleflux);
/******************************** buildmultisky.c ****************************/
/* FITS reading, multiple files
@ -299,9 +300,10 @@ write_world_coords(const char *imgfile, GHashTable *pixtable, double minpix, dou
bmaj (rad),bmin (rad) ,pa (rad): PSF parameters from the FITS files: Nfx1 arrays
minpix: footprint of MEAN psf in pixels
ignlist: list of islands to ignore (integer list)
donegative: if >0, fit -ve pixels instead of +ve pixels
*/
extern int
read_fits_file_f(const char *fitsdir, const char *maskfile, GHashTable **pixtable, int *Nf, double **freqs, double **bmaj, double **bmin, double **pa, int beam_given, double *minpix, GList *ignlist);
read_fits_file_f(const char *fitsdir, const char *maskfile, GHashTable **pixtable, int *Nf, double **freqs, double **bmaj, double **bmin, double **pa, int beam_given, double *minpix, GList *ignlist, int donegative);
/* calculates the centroids, etc
@ -353,9 +355,11 @@ filter_pixels_f(GHashTable *pixtable, double wcutoff);
clusterratio: components closer than clusteratio*(bmaj+bmin)/2 will be merged
nclusters: no of max clusters to cluster the sky
unistr: unique string for source names
donegative: if >0, print -ve flux instead of +ve flux
scaleflux: if 1, scale model flux to match total flux of island
*/
extern int
write_world_coords_f(const char *imgfile, GHashTable *pixtable, double minpix, int Nf, double *freqs, double *bmaj, double *bmin, double *bpa, double ref_freq, int outformat, double clusterratio, int nclusters, const char *unistr);
write_world_coords_f(const char *imgfile, GHashTable *pixtable, double minpix, int Nf, double *freqs, double *bmaj, double *bmin, double *bpa, double ref_freq, int outformat, double clusterratio, int nclusters, const char *unistr,int donegative, int scaleflux);
/******************************** fitpixels.c *******************************/

148
src/buildsky/convert_skymodel.py
View File

@ -62,7 +62,7 @@ def convert_sky_bbs_lsm(infilename,outfilename):
\s* # skip white space
\, # skip comma
\s* # skip white space
(?P<col10>[-+]?\d+(\.\d+)?) # sI
(?P<col10>[-+]?\d+(\.\d+)?([eE](-|\+)(\d+))?) # sI
\s* # skip white space
\, # skip comma
\s* # skip white space
@ -106,7 +106,7 @@ def convert_sky_bbs_lsm(infilename,outfilename):
\s* # skip white space
\, # skip comma
\s* # skip white space
(?P<col10>[-+]?\d+(\.\d+)?) # sI
(?P<col10>[-+]?\d+(\.\d+)?([eE](-|\+)(\d+))?) # sI
\s* # skip white space
\, # skip comma
\s* # skip white space
@ -274,54 +274,9 @@ def convert_sky_bbs_lsm(infilename,outfilename):
outfile.write("## LSM file\n")
outfile.write("### Name | RA (hr,min,sec) | DEC (deg,min,sec) | I | Q | U | V | SI | RM | eX | eY | eP | freq0\n")
for eachline in all:
v=pp1.search(eachline)
v=pp3.search(eachline)
if v!= None:
strline=str(v.group('col1'))+' '+str(v.group('col4'))+' '+str(v.group('col5'))+' '+str(v.group('col6'))
strline=strline+' '+str(v.group('col7'))+' '+str(v.group('col8'))+' '+str(v.group('col9'))+' '+str(v.group('col10'))
strline=strline+' '+str(v.group('col11'))+' '+str(v.group('col12'))+' '+str(v.group('col13'))+' '+str(v.group('col15'))
strline=strline+' 0 0 0 0 '+str(v.group('col14'))+'\n'
outfile.write(strline)
else:
v=pp.search(eachline)
if v!= None:
strline=str(v.group('col1'))+' '+str(v.group('col4'))+' '+str(v.group('col5'))+' '+str(v.group('col6'))
strline=strline+' '+str(v.group('col7'))+' '+str(v.group('col8'))+' '+str(v.group('col9'))+' '+str(v.group('col10'))
strline=strline+' 0 0 0 0 0 0 0 0\n'
outfile.write(strline)
else:
v=pp2.search(eachline)
if v!= None:
stype=v.group('col2')
sname=str(v.group('col1'))
bad_source=False
if stype=='GAUSSIAN':
sname='G'+sname
strline=sname+' '+str(v.group('col4'))+' '+str(v.group('col5'))+' '+str(v.group('col6'))
strline=strline+' '+str(v.group('col7'))+' '+str(v.group('col8'))+' '+str(v.group('col9'))+' '+str(v.group('col10'))
strline=strline+' '+str(v.group('col11'))+' '+str(v.group('col12'))+' '+str(v.group('col13'))+' '+str(v.group('col18'))
# need the right conversion factor for Gaussians: Bmaj,Bmin arcsec (rad) diameter, PA (deg) West-> counterclockwise, BDSM its North->counterclockwise
if stype=='GAUSSIAN':
bmaj=float(v.group('col14'))*(0.5/3600)*math.pi/180.0
bmin=float(v.group('col15'))*(0.5/3600)*math.pi/180.0
bpa=math.pi/2-(math.pi-float(v.group('col16'))/180.0*math.pi)
# also throw away bad Gaussians with zero bmaj or bmin
if bmaj<1e-6 or bmin<1e-6:
bad_source=True
else:
bmaj=0.0
bmin=0.0
bpa=0.0
strline=strline+' 0 '+str(bmaj)+' '+str(bmin)+' '+str(bpa)
strline=strline+' '+str(v.group('col17'))+'\n'
# only write good sources
if not bad_source:
outfile.write(strline)
else:
print 'Error in source '+strline
else:
v=pp3.search(eachline)
if v!= None:
################################################################################
stype=v.group('col2')
sname=str(v.group('col1'))
bad_source=False
@ -349,35 +304,72 @@ def convert_sky_bbs_lsm(infilename,outfilename):
if not bad_source:
outfile.write(strline)
else:
print 'Error in source '+strline
else:
splitfields=eachline.split(',');
if len(splitfields)>10:
# type
if 'GAUSSIAN' in splitfields[1]:
# read last 3 fields for shape
bpa=splitfields.pop()
bmin=splitfields.pop()
bmaj=splitfields.pop()
# convert to right format
bmaj=float(bmaj)*(0.5/3600.0)*math.pi/180.0
bmin=float(bmin)*(0.5/3600.0)*math.pi/180.0
bpa=math.pi/2-(math.pi-float(bpa)/180.0*math.pi)
sname='G'+splitfields[0]
else:
sname='P'+splitfields[0]
bmaj=0
bmin=0
bpa=0
# read in 2nd,3rd field for RA,DEC
raline=splitfields[2].split(':')
decline=splitfields[3].split('.')
# finally sI
sI=splitfields[4]
strline=sname+' '+str(int(raline[0]))+' '+str(int(raline[1]))+' '+str(float(raline[2]))+' '+str(int(decline[0]))+' '+str(int(decline[1]))+' '+str(float(decline[2]+'.'+decline[3]))+' '+str(float(sI))
# add missing fields (Q U V SI RM)
strline=strline+' 0 0 0 0 0 '+str(bmaj)+' '+str(bmin)+' '+str(bpa)+' 150e6\n'
outfile.write(strline)
pass
#print 'Error in source '+strline
################################################################################
else:
v=pp2.search(eachline)
if v!= None:
################################################################################
print eachline
stype=v.group('col2')
sname=str(v.group('col1'))
bad_source=False
if stype=='GAUSSIAN':
sname='G'+sname
strline=sname+' '+str(v.group('col4'))+' '+str(v.group('col5'))+' '+str(v.group('col6'))
strline=strline+' '+str(v.group('col7'))+' '+str(v.group('col8'))+' '+str(v.group('col9'))+' '+str(v.group('col10'))
strline=strline+' '+str(v.group('col11'))+' '+str(v.group('col12'))+' '+str(v.group('col13'))+' '+str(v.group('col18'))
# need the right conversion factor for Gaussians: Bmaj,Bmin arcsec (rad) diameter, PA (deg) West-> counterclockwise, BDSM its North->counterclockwise
if stype=='GAUSSIAN':
bmaj=float(v.group('col14'))*(0.5/3600)*math.pi/180.0
bmin=float(v.group('col15'))*(0.5/3600)*math.pi/180.0
bpa=math.pi/2-(math.pi-float(v.group('col16'))/180.0*math.pi)
# also throw away bad Gaussians with zero bmaj or bmin
if bmaj<1e-6 or bmin<1e-6:
bad_source=True
print "%f %f %f"%(bmaj,bmin,bpa)
else:
bmaj=0.0
bmin=0.0
bpa=0.0
strline=strline+' 0 '+str(bmaj)+' '+str(bmin)+' '+str(bpa)
strline=strline+' '+str(v.group('col17'))+'\n'
print strline
# only write good sources
if not bad_source:
outfile.write(strline)
else:
pass
#print 'Error in source '+strline
################################################################################
else:
v=pp1.search(eachline)
if v!= None:
################################################################################
strline=str(v.group('col1'))+' '+str(v.group('col4'))+' '+str(v.group('col5'))+' '+str(v.group('col6'))
strline=strline+' '+str(v.group('col7'))+' '+str(v.group('col8'))+' '+str(v.group('col9'))+' '+str(v.group('col10'))
strline=strline+' '+str(v.group('col11'))+' '+str(v.group('col12'))+' '+str(v.group('col13'))+' '+str(v.group('col15'))
strline=strline+' 0 0 0 0 '+str(v.group('col14'))+'\n'
outfile.write(strline)
################################################################################
else:
v=pp.search(eachline)
if v!= None:
################################################################################
strline=str(v.group('col1'))+' '+str(v.group('col4'))+' '+str(v.group('col5'))+' '+str(v.group('col6'))
strline=strline+' '+str(v.group('col7'))+' '+str(v.group('col8'))+' '+str(v.group('col9'))+' '+str(v.group('col10'))
strline=strline+' 0 0 0 0 0 0 0 0\n'
outfile.write(strline)
################################################################################
else:
pass
#print eachline
def convert_sky_lsm_bbs(infilename,outfilename):
# LSM format

29
src/buildsky/create_clusters.py
View File

@ -115,15 +115,34 @@ def read_lsm_sky(infilename):
for eachline in all:
v=pp1.search(eachline)
if v!= None:
# find RA,DEC (rad) and flux
mra=(float(v.group('col2'))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0*math.pi/180.0
mdec=(float(v.group('col5'))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)*math.pi/180.0
# find RA,DEC (rad) and flux, with proper sign
if (float(v.group('col2'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mra=mysign*(abs(float(v.group('col2')))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0*math.pi/180.0
if (float(v.group('col5'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mdec=mysign*(abs(float(v.group('col5')))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)*math.pi/180.0
SR[str(v.group('col1'))]=(v.group('col1'),mra,mdec,float(v.group('col8')))
else:
v=pp.search(eachline)
if v!= None:
mra=(float(v.group('col2'))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0*math.pi/180.0
mdec=(float(v.group('col5'))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)*math.pi/180.0
if (float(v.group('col2'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mra=mysign*(abs(float(v.group('col2')))+float(v.group('col3'))/60.0+float(v.group('col4'))/3600.0)*360.0/24.0*math.pi/180.0
if (float(v.group('col5'))) >= 0.0:
mysign=1.0
else:
mysign=-1.0
mdec=mysign*(abs(float(v.group('col5')))+float(v.group('col6'))/60.0+float(v.group('col7'))/3600.0)*math.pi/180.0
SR[str(v.group('col1'))]=(v.group('col1'),mra,mdec,float(v.group('col8')))
print 'Read %d sources'%len(SR)

18
src/buildsky/fitmultipixels.c
View File

@ -266,7 +266,7 @@ mylm_fit_single_pf_1d(double *p, double *x, int m, int n, void *data) {
double
fit_single_point_f(hpixelf *parr, int npix, int Nf, double *freqs, double *bmaj, double *bmin, double *bpa, double ref_freq, int maxiter, double *ll1, double *mm1, double *sI1, double *sP1){
int ci,cj,ret;
int ci,cj;
double *p,*p1,*p2, // params m x 1
*x; // observed data n x 1, the image pixel fluxes
int m,n;
@ -381,12 +381,12 @@ fit_single_point_f(hpixelf *parr, int npix, int Nf, double *freqs, double *bmaj,
double aic1,aic2,aic3;
//ret=dlevmar_dif(mylm_fit_single_sipf, p, x, m, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
ret=clevmar_der_single_nocuda(mylm_fit_single_sipf, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single_sipf, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
/* penalize only 1/10 of parameters */
aic3=0.3+log(info[1]);
ret=clevmar_der_single_nocuda(mylm_fit_single_sipf_2d, NULL, p2, x, m-1, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single_sipf_2d, NULL, p2, x, m-1, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
aic2=0.2+log(info[1]);
ret=clevmar_der_single_nocuda(mylm_fit_single_sipf_1d, NULL, p1, x, m-2, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single_sipf_1d, NULL, p1, x, m-2, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
aic1=0.1+log(info[1]);
/* choose one with minimum error */
if (aic3<aic2) {
@ -467,6 +467,8 @@ fit_single_point_f(hpixelf *parr, int npix, int Nf, double *freqs, double *bmaj,
sP1[2]=p[3];
free(p);
free(p1);
free(p2);
free(x);
/* AIC, 4 parms */
@ -479,7 +481,7 @@ fit_single_point_f(hpixelf *parr, int npix, int Nf, double *freqs, double *bmaj,
double
fit_N_point_em_f(hpixelf *parr, int npix, int Nf, double *freqs, double *bmaj, double *bmin, double *bpa, double ref_freq, int maxiter, int max_em_iter, double *ll, double *mm, double *sI, double *sP, int N, int Nh, hpoint *hull){
int ci,ret,cj,ck;
int ci,cj,ck;
double *p,*p1,*p2, // params m x 1
*x; // observed data n x 1, the image pixel fluxes
double *xdummy, *xsub; //extra arrays
@ -647,12 +649,12 @@ fit_N_point_em_f(hpixelf *parr, int npix, int Nf, double *freqs, double *bmaj, d
p[3]=sP[cj+2*N]; p1[3]=p2[3]=0.0;
//ret=dlevmar_dif(mylm_fit_single_pf, p, xdummy, m, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
ret=clevmar_der_single_nocuda(mylm_fit_single_pf, NULL, p, xdummy, m, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single_pf, NULL, p, xdummy, m, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
/* penalize only 1/10 of parameters */
aic3=0.3+log(info[1]);
ret=clevmar_der_single_nocuda(mylm_fit_single_pf_2d, NULL, p2, xdummy, m-1, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single_pf_2d, NULL, p2, xdummy, m-1, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
aic2=0.2+log(info[1]);
ret=clevmar_der_single_nocuda(mylm_fit_single_pf_1d, NULL, p1, xdummy, m-2, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single_pf_1d, NULL, p1, xdummy, m-2, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
aic1=0.1+log(info[1]);
/* choose one with minimum error */
if (aic3<aic2) {

18
src/buildsky/fitpixels.c
View File

@ -293,7 +293,7 @@ mylm_fit_single(double *p, double *x, int m, int n, void *data) {
double
fit_single_point(hpixel *parr, int npix, double bmaj, double bmin, double bpa, int maxiter, double *ll1, double *mm1, double *sI1){
int ci,ret;
int ci;
double *p, // params m x 1
*x; // observed data n x 1, the image pixel fluxes
int m,n;
@ -343,8 +343,8 @@ fit_single_point(hpixel *parr, int npix, double bmaj, double bmin, double bpa, i
lmdata.bmin=bmin;
lmdata.bpa=bpa;
//ret=dlevmar_dif(mylm_fit_single, p, x, 3, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
ret=clevmar_der_single_nocuda(mylm_fit_single, NULL, p, x, 3, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
//dlevmar_dif(mylm_fit_single, p, x, 3, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_single, NULL, p, x, 3, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
#ifdef DEBUG
print_levmar_info(info[0],info[1],(int)info[5], (int)info[6], (int)info[7], (int)info[8], (int)info[9]);
@ -398,7 +398,7 @@ mylm_fit_N(double *p, double *x, int m, int n, void *data) {
double
fit_N_point_em(hpixel *parr, int npix, double bmaj, double bmin, double bpa, int maxiter, int max_em_iter, double *ll, double *mm, double *sI, int N, int Nh, hpoint *hull){
int ci,ret,cj,ck;
int ci,cj,ck;
double *p, // params m x 1
*x; // observed data n x 1, the image pixel fluxes
double *xdummy, *xsub; //extra arrays
@ -493,17 +493,17 @@ fit_N_point_em(hpixel *parr, int npix, double bmaj, double bmin, double bpa, int
my_daxpy(n, xsub, -1.0, xdummy);
}
}
//ret=dlevmar_dif(mylm_fit_single, &p[3*cj], xdummy, 3, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
//ret=clevmar_der_single_nocuda(mylm_fit_single, NULL, &p[3*cj], xdummy, 3, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
//dlevmar_dif(mylm_fit_single, &p[3*cj], xdummy, 3, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
//clevmar_der_single_nocuda(mylm_fit_single, NULL, &p[3*cj], xdummy, 3, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
// weighted least squares estimation
ret=dweighted_ls_fit(&p[3*cj], xdummy, 3, n, maxiter, (void*)&lmdata);
dweighted_ls_fit(&p[3*cj], xdummy, 3, n, maxiter, (void*)&lmdata);
}
}
/* final global fit */
//ret=dlevmar_dif(mylm_fit_N, p, x, m, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
ret=clevmar_der_single_nocuda(mylm_fit_N, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
//dlevmar_dif(mylm_fit_N, p, x, m, n, maxiter, opts, info, NULL, NULL, (void*)&lmdata); // no Jacobian
clevmar_der_single_nocuda(mylm_fit_N, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata); // no Jacobian
#ifdef DEBUG
print_levmar_info(info[0],info[1],(int)info[5], (int)info[6], (int)info[7], (int)info[8], (int)info[9]);
printf("Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);

26
src/buildsky/main.c
View File

@ -34,7 +34,7 @@ print_help(void) {
fprintf(stderr,"PSF options\n-a : major axis width (arcsec): default 10\n");
fprintf(stderr,"-b : minor axis width (arcsec): default 10\n");
fprintf(stderr,"-p : position angle, measured from positive y axis, counter clockwise (deg): default 0\n");
fprintf(stderr,"If NO PSF is given, it will be read from the FITS file\n");
fprintf(stderr,"If NO PSF is given, it will be read from the FITS file\n\n");
fprintf(stderr,"-o format: output format (0: BBS, 1: LSM (with 3 order spec.idx) 2: (l,m) positions) default BBS\n");
fprintf(stderr,"-g ignorelist: text file of island numbers to ignore: default none\n");
fprintf(stderr,"-w cutoff: cutoff value to detect possible sidelobes: defalut 0\n");
@ -42,13 +42,15 @@ print_help(void) {
fprintf(stderr,"-l maxfits: limit the number of attempted fits to this value: default 10\n");
fprintf(stderr,"-k clusters: max no of clusters to cluster the sky\n(-ve clusters for hierarchical, 0 for no clustering)\n");
fprintf(stderr,"-s string: additional unique string to use with source names: default nothing\n");
fprintf(stderr,"-N : if used, fit negative flux instead of positive\n");
fprintf(stderr,"-q : if 1, scale model fluxes to match the total flux of detected island: default 0\n");
fprintf(stderr,"Report bugs to <sarod@users.sf.net>\n");
}
void
print_copyright(void) {
printf("Buildsky 0.0.9 (C) 2011-2016 Sarod Yatawatta\n");
printf("Buildsky 0.1.0 (C) 2011-2016 Sarod Yatawatta\n");
}
@ -68,6 +70,7 @@ int main(int argc, char **argv) {
double threshold=0.0;
int use_em=1;
int donegative=0;
int maxiter=1000;
int maxemiter=4;
int outformat=0; /* 0: BBS, 1: LSM */
@ -85,6 +88,8 @@ int main(int argc, char **argv) {
int nclusters=0;
double wcutoff=0.0;
int scaleflux=0; /* if 1, scale flux to match the total flux of island */
/* for multiple FITS files */
int Nf=0;
double *freqs,*bmajs,*bmins,*bpas;
@ -97,7 +102,7 @@ int main(int argc, char **argv) {
print_help();
return 1;
}
while ((c=getopt(argc,argv,"f:m:t:i:e:a:b:o:p:l:g:c:s:d:k:w:nh"))!=-1) {
while ((c=getopt(argc,argv,"a:b:c:d:e:f:g:i:k:l:m:o:p:q:s:t:w:Nnh"))!=-1) {
switch(c) {
case 'f':
if (optarg) {
@ -190,6 +195,12 @@ int main(int argc, char **argv) {
if (endptr==optarg) maxemiter=1000;
}
break;
case 'q':
if (optarg) {
scaleflux=(int)strtol(optarg,&endptr,base);
if (endptr==optarg) scaleflux=0;
}
break;
case 'a':
if (optarg) {
bmaj=strtod(optarg,0);
@ -226,6 +237,9 @@ int main(int argc, char **argv) {
case 'n':
use_em=0;
break;
case 'N':
donegative=1;
break;
case 'w':
if (optarg) {
wcutoff=strtod(optarg,0);
@ -279,7 +293,7 @@ int main(int argc, char **argv) {
read_fits_file(ffile,mfile, &pixtable,&bmaj,&bmin,&bpa, beam_given, &minpix, ignlist);
} else {
freqs=bmajs=bmins=bpas=0;
read_fits_file_f(ffile,mfile, &pixtable,&Nf,&freqs,&bmajs,&bmins,&bpas, beam_given, &minpix, ignlist);
read_fits_file_f(ffile,mfile, &pixtable,&Nf,&freqs,&bmajs,&bmins,&bpas, beam_given, &minpix, ignlist,donegative);
}
/* dont need ignore list anymore */
for(ignli=ignlist; ignli!=NULL; ignli=g_list_next(ignli)) {
@ -300,10 +314,10 @@ int main(int argc, char **argv) {
}
if (!multifits) {
process_pixels(pixtable,bmaj,bmin,bpa,minpix,threshold,maxiter,maxemiter,use_em,maxfits);
write_world_coords(ffile,pixtable,minpix,bmaj,bmin,bpa,outformat,clusterratio,nclusters,unistr);
write_world_coords(ffile,pixtable,minpix,bmaj,bmin,bpa,outformat,clusterratio,nclusters,unistr,scaleflux);
} else {
process_pixels_f(pixtable,Nf,freqs,bmajs,bmins,bpas,&ref_freq,minpix,threshold,maxiter,maxemiter,use_em,maxfits);
write_world_coords_f(mfile,pixtable,minpix,Nf,freqs,bmajs,bmins,bpas,ref_freq,outformat,clusterratio,nclusters,unistr);
write_world_coords_f(mfile,pixtable,minpix,Nf,freqs,bmajs,bmins,bpas,ref_freq,outformat,clusterratio,nclusters,unistr,donegative, scaleflux);
free(freqs);
free(bmajs);
free(bmins);

8
src/buildsky/scluster.c
View File

@ -99,7 +99,7 @@ fit_two_to_one(GList *pixlist, double bmaj, double bmin, double bpa, double *lva
hpixel *ppix;
double *pixell,*pixelm,*p,*x,*x1;
GList *pli;
int ci,ret;
int ci;
/* setup levmar */
double opts[CLM_OPTS_SZ], info[CLM_INFO_SZ];
opts[0]=CLM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-15;
@ -177,7 +177,7 @@ fit_two_to_one(GList *pixlist, double bmaj, double bmin, double bpa, double *lva
printf("Initial merge val (%lf)\n",p[0]);
#endif
ret=clevmar_der_single_nocuda(mylm_fit_single_pfmult, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata);
clevmar_der_single_nocuda(mylm_fit_single_pfmult, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata);
/* output */
*sI=p[0];
@ -252,7 +252,7 @@ fit_two_to_one_f(GList *pixlist, int Nf, double *freqs, double *bmaj, double *bm
hpixelf *ppix;
double *pixell,*pixelm,*p,*x,*x1;
GList *pli;
int ci,ret;
int ci;
/* setup levmar */
double opts[CLM_OPTS_SZ], info[CLM_INFO_SZ];
opts[0]=CLM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-15;
@ -338,7 +338,7 @@ fit_two_to_one_f(GList *pixlist, int Nf, double *freqs, double *bmaj, double *bm
printf("Initial merge val (%lf,%lf)\n",p[0],p[1]);
#endif
ret=clevmar_der_single_nocuda(mylm_fit_single_sipfmult, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata);
clevmar_der_single_nocuda(mylm_fit_single_sipfmult, NULL, p, x, m, n, maxiter, opts, info, 2, (void*)&lmdata);
/* output */
*sI=p[0];