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

update

This commit is contained in:
Sarod Yatawatta 2018-02-07 10:55:51 +01:00
parent 3a95d1012b
commit 295e6f9c4b
3 changed files with 178 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
src/MS/data.cpp
View File

@ -18,7 +18,6 @@
*/
#include "data.h"
// #include "sagecal.h"
#include "Radio.h"
#include <measures/Measures/MDirection.h>
#include <measures/Measures/UVWMachine.h>
@ -87,6 +86,13 @@ int Data::Nskip=0;
int Data::verbose=0; /* no verbose output */
int Data::mdl=0; /* no AIC/MDL calculation by default */
int Data::GPUpredict=0; /* use CPU for model calculation, if GPU not specified */
#ifdef HAVE_CUDA
int Data::heapsize=GPU_HEAP_SIZE; /* heap size in GPU (MB) to be used in malloc() */
#endif
int Data::servermode=-1; /* by default, no client-server mode */
char *Data::servername=NULL;
char *Data::portnumber=NULL;
using namespace Data;
@ -141,6 +147,7 @@ Data::readAuxData(const char *fname, Data::IOData *data) {
data->Nchan=chan_freq.shape(0)[0];
data->Nms=1;
/* allocate memory */
try {
data->u=new double[data->Nbase*data->tilesz];
data->v=new double[data->Nbase*data->tilesz];
data->w=new double[data->Nbase*data->tilesz];
@ -149,6 +156,10 @@ Data::readAuxData(const char *fname, Data::IOData *data) {
data->freqs=new double[data->Nchan];
data->flag=new double[data->Nbase*data->tilesz];
data->NchanMS=new int[data->Nms];
} catch (const std::bad_alloc& e) {
cout<<"Allocating memory for data failed. Quitting."<< e.what() << endl;
exit(1);
}
data->NchanMS[0]=data->Nchan;
/* copy freq */
@ -192,6 +203,7 @@ Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
data->Nchan=chan_freq.shape(0)[0];
data->Nms=1;
/* allocate memory */
try {
data->u=new double[data->Nbase*data->tilesz];
data->v=new double[data->Nbase*data->tilesz];
data->w=new double[data->Nbase*data->tilesz];
@ -200,6 +212,10 @@ Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
data->freqs=new double[data->Nchan];
data->flag=new double[data->Nbase*data->tilesz];
data->NchanMS=new int[data->Nms];
} catch (const std::bad_alloc& e) {
cout<<"Allocating memory for data failed. Quitting."<< e.what() << endl;
exit(1);
}
data->NchanMS[0]=data->Nchan;
/* copy freq */
@ -522,7 +538,7 @@ Data::readAuxDataList(vector<string> msnames, Data::IOData *data) {
/* each time this is called read in data from MS, and format them as
u,v,w: u,v,w coordinates (wavelengths) size Nbase*tilesz x 1
u,v,w are ordered with baselines, timeslots
x: data to write size Nbase*8*tileze x 1
x: data to write size Nbase*8*tilesz x 1
ordered by XX(re,im),XY(re,im),YX(re,im), YY(re,im), baseline, timeslots
fratio: flagged data as a ratio to all available data
*/
@ -544,8 +560,8 @@ Data::loadData(Table ti, Data::IOData iodata, double *fratio) {
/* check we get correct rows */
int nrow=t.nrow();
if(nrow-iodata.N*iodata.tilesz>iodata.tilesz*iodata.Nbase) {
cout<<"Error in rows"<<endl;
if(nrow<iodata.tilesz*iodata.Nbase-iodata.tilesz*iodata.N) {
cout<<"Warning: Missing rows, got "<<nrow<<" expect "<<iodata.tilesz*iodata.Nbase<<" +- "<<iodata.tilesz*iodata.N<<". (probably the last time interval, so not a big issue)."<<endl;
}
int row0=0;
/* tapering */
@ -558,7 +574,7 @@ Data::loadData(Table ti, Data::IOData iodata, double *fratio) {
}
/* counters for finding flagged data ratio */
int countgood=0; int countbad=0;
for(int row = 0; row < nrow; row++) {
for(int row = 0; row < nrow && row0<iodata.tilesz*iodata.Nbase; row++) {
uInt i = a1(row); //antenna1
uInt j = a2(row); //antenna2
/* only work with cross correlations */
@ -651,6 +667,16 @@ Data::loadData(Table ti, Data::IOData iodata, double *fratio) {
if (row0<iodata.tilesz*iodata.Nbase) {
for(int row = row0; row<iodata.tilesz*iodata.Nbase; row++) {
iodata.flag[row]=1;
}
/* set uvw and data to 0 to eliminate any funny business */
memset(&iodata.u[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
memset(&iodata.v[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
memset(&iodata.w[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
memset(&iodata.x[8*row0],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
for(int k = 0; k < iodata.Nchan; k++) {
memset(&iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
}
}
/* flagged data / total usable data, not counting excluded baselines */
@ -689,8 +715,8 @@ Data::loadData(Table ti, Data::IOData iodata, LBeam binfo, double *fratio) {
/* check we get correct rows */
int nrow=t.nrow();
if(nrow-iodata.N*iodata.tilesz>iodata.tilesz*iodata.Nbase) {
cout<<"Error in rows"<<endl;
if(nrow<iodata.tilesz*iodata.Nbase-iodata.tilesz*iodata.N) {
cout<<"Warning: Missing rows, got "<<nrow<<" expect "<<iodata.tilesz*iodata.Nbase<<" +- "<<iodata.tilesz*iodata.N<<". (probably the last time interval, so not a big issue)."<<endl;
}
int row0=0;
int rowt=0;
@ -704,13 +730,11 @@ Data::loadData(Table ti, Data::IOData iodata, LBeam binfo, double *fratio) {
}
/* counters for finding flagged data ratio */
int countgood=0; int countbad=0;
for(int row = 0; row < nrow; row++) {
for(int row = 0; row < nrow && row0<iodata.tilesz*iodata.Nbase; row++) {
uInt i = a1(row); //antenna1
uInt j = a2(row); //antenna2
if (!i && !j) {/* use baseline 0-0 to extract time */
double tt=tut(row);
//cout.precision(22);
//cout<<"("<<i<<","<<j<<") "<<rowt<<"="<<tt<<endl;
/* convert MJD (s) to JD (days) */
binfo.time_utc[rowt++]=(tt/86400.0+2400000.5); /* no +0.5 added */
}
@ -804,6 +828,16 @@ Data::loadData(Table ti, Data::IOData iodata, LBeam binfo, double *fratio) {
if (row0<iodata.tilesz*iodata.Nbase) {
for(int row = row0; row<iodata.tilesz*iodata.Nbase; row++) {
iodata.flag[row]=1;
}
/* set uvw and data to 0 to eliminate any funny business */
memset(&iodata.u[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
memset(&iodata.v[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
memset(&iodata.w[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
memset(&iodata.x[8*row0],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
for(int k = 0; k < iodata.Nchan; k++) {
memset(&iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
}
}
/* flagged data / total usable data, not counting excluded baselines */
@ -993,7 +1027,7 @@ Data::writeData(Table ti, Data::IOData iodata) {
/* check we get correct rows */
int nrow=t.nrow();
if(nrow-iodata.N*iodata.tilesz>iodata.tilesz*iodata.Nbase) {
cout<<"Error in rows"<<endl;
cout<<"Warning: Missing rows, got "<<nrow<<" expect "<<iodata.tilesz*iodata.Nbase<<" +- "<<iodata.tilesz*iodata.N<<". (probably the last time interval, so not a big issue)."<<endl;
}
//cout<<"Table rows "<<nrow<<" Data rows "<<iodata.tilesz*iodata.Nbase+iodata.tilesz*iodata.N<<endl;
int row0=0;

5
src/MS/data.h
View File

@ -165,5 +165,10 @@ namespace Data
extern int verbose; /* if >0, enable verbose output */
extern int mdl; /* if given, calculate AIC/MDL for different poly configs and find minimum */
extern int GPUpredict; /* if given, use GPU for model calculation */
extern int heapsize; /* heap size in GPU (MB), for using malloc() */
/* for client server mode */
extern int servermode; /* 0: client, 1: server, else default operation */
extern char *servername; /* server host name or ip address */
extern char *portnumber; /* which port number to use for communication */
}
#endif //__DATA_H__

144
src/MS/main.cpp
View File

@ -25,8 +25,6 @@
#include <pthread.h>
#include <casacore/casa/Quanta/Quantum.h>
#include "cuda_profiler_api.h"
#include <Dirac.h>
#include <Radio.h>
@ -68,7 +66,9 @@ print_help(void) {
cout << "-z ignore_clusters: if only doing a simulation, ignore the cluster ids listed in this file" << endl;
cout << "-b 0,1 : if 1, solve for each channel: default " <<Data::doChan<< endl;
cout << "-B 0,1 : if 1, predict array beam: default " <<Data::doBeam<< endl;
#ifdef HAVE_CUDA
cout << "-E 0,1 : if 1, use GPU for model computing: default " <<Data::GPUpredict<< endl;
#endif
cout << "-x exclude baselines length (lambda) lower than this in calibration : default "<<Data::min_uvcut << endl;
cout << "-y exclude baselines length (lambda) higher than this in calibration : default "<<Data::max_uvcut << endl;
cout <<endl<<"Advanced options:"<<endl;
@ -80,6 +80,9 @@ print_help(void) {
cout << "-H robust nu, upper bound: default "<<Data::nuhigh<< endl;
cout << "-W pre-whiten data: default "<<Data::whiten<< endl;
cout << "-R randomize iterations: default "<<Data::randomize<< endl;
#ifdef HAVE_CUDA
cout << "-S GPU heap size (MB): default "<<Data::heapsize<< endl;
#endif
cout << "-D 0,1,2 : if >0, enable diagnostics (Jacobian Leverage) 1 replace Jacobian Leverage as output, 2 only fractional noise/leverage is printed: default " <<Data::DoDiag<< endl;
cout << "-q solutions.txt: if given, initialize solutions by reading this file (need to have the same format as a solution file, only solutions for 1 timeslot needed)"<< endl;
cout <<"Report bugs to <sarod@users.sf.net>"<<endl;
@ -95,7 +98,7 @@ ParseCmdLine(int ac, char **av) {
print_help();
exit(0);
}
while((c=getopt(ac, av, "a:b:c:d:e:f:g:j:k:l:m:n:o:p:q:s:t:x:y:z:B:D:E:F:I:J:O:L:H:R:W:E:h"))!= -1)
while((c=getopt(ac, av, "a:b:c:d:e:f:g:j:k:l:m:n:o:p:q:s:t:x:y:z:B:D:E:F:I:J:O:L:H:R:S:W:E:h"))!= -1)
{
switch(c)
{
@ -184,6 +187,11 @@ ParseCmdLine(int ac, char **av) {
case 'J':
phaseOnly= atoi(optarg);
break;
#ifdef HAVE_CUDA
case 'S':
heapsize= atoi(optarg);
break;
#endif
case 'x':
Data::min_uvcut= atof(optarg);
break;
@ -244,7 +252,11 @@ main(int argc, char **argv) {
Data::readMSlist(Data::MSlist,&msnames);
}
if (Data::TableName) {
if (!doBeam) {
Data::readAuxData(Data::TableName,&iodata);
} else {
Data::readAuxData(Data::TableName,&iodata,&beam);
}
cout<<"Only one MS"<<endl;
} else if (Data::MSlist) {
Data::readAuxDataList(msnames,&iodata);
@ -255,12 +267,11 @@ main(int argc, char **argv) {
srand(time(0)); /* use different seed */
}
// openblas_set_num_threads(1);//Data::Nt;
// export OMP_NUM_THREADS=1
openblas_set_num_threads(1);//Data::Nt;
/**********************************************************/
int M,Mt,ci,cj,ck;
/* parameters */
double *p,*pinit;
double *p,*pinit,*pfreq;
double **pm;
complex double *coh;
FILE *sfp=0;
@ -333,6 +344,19 @@ main(int argc, char **argv) {
}
}
#ifdef USE_MIC
/* need for bitwise copyable parameter passing */
int *mic_pindex,*mic_chunks;
if ((mic_chunks=(int*)calloc((size_t)M,sizeof(int)))==0) {
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
exit(1);
}
if ((mic_pindex=(int*)calloc((size_t)Mt,sizeof(int)))==0) {
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
exit(1);
}
int cl=0;
#endif
/* update cluster array with correct pointers to parameters */
cj=0;
@ -341,8 +365,14 @@ main(int argc, char **argv) {
fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
exit(1);
}
#ifdef USE_MIC
mic_chunks[ci]=carr[ci].nchunk;
#endif
for (ck=0; ck<carr[ci].nchunk; ck++) {
carr[ci].p[ck]=cj*8*iodata.N;
#ifdef USE_MIC
mic_pindex[cl++]=carr[ci].p[ck];
#endif
cj++;
}
}
@ -399,8 +429,8 @@ main(int argc, char **argv) {
double res_0,res_1,res_00,res_01;
/* previous residual */
// double res_prev=CLM_DBL_MAX;
// double res_ratio=5; /* how much can the residual increase before resetting solutions */
double res_prev=CLM_DBL_MAX;
double res_ratio=5; /* how much can the residual increase before resetting solutions */
res_0=res_1=res_00=res_01=0.0;
/**********************************************************/
@ -443,18 +473,29 @@ main(int argc, char **argv) {
/* starting iterations are doubled */
// int start_iter=1;
// int sources_precessed=0;
int start_iter=1;
int sources_precessed=0;
double inv_c=1.0/CONST_C;
#ifdef HAVE_CUDA
cudaProfilerStart();
/* setup Heap of GPU, only need to be done once, before any kernel is launched */
if (GPUpredict>0) {
for (int gpuid=0; gpuid<=MAX_GPU_ID; gpuid++) {
cudaSetDevice(gpuid);
cudaDeviceSetLimit(cudaLimitMallocHeapSize, Data::heapsize*1024*1024);
}
}
#endif
while (msitr[0]->more()) {
start_time = time(0);
if (iodata.Nms==1) {
if (!doBeam) {
Data::loadData(msitr[0]->table(),iodata,&iodata.fratio);
} else {
Data::loadData(msitr[0]->table(),iodata,beam,&iodata.fratio);
}
} else {
Data::loadDataList(msitr,iodata,&iodata.fratio);
}
@ -469,11 +510,41 @@ main(int argc, char **argv) {
preset_flags_and_data(iodata.Nbase*iodata.tilesz,iodata.flag,barr,iodata.x,Data::Nt);
/* if data is being whitened, whiten x here,
no need for a copy because we use xo for residual calculation */
if (Data::whiten) {
whiten_data(iodata.Nbase*iodata.tilesz,iodata.x,iodata.u,iodata.v,iodata.freq0,Data::Nt);
}
/* precess source locations (also beam pointing) from J2000 to JAPP if we do any beam predictions,
using first time slot as epoch */
// sources_precessed=1;
if (doBeam && !sources_precessed) {
precess_source_locations(beam.time_utc[iodata.tilesz/2],carr,M,&beam.p_ra0,&beam.p_dec0,Data::Nt);
sources_precessed=1;
}
#ifdef USE_MIC
double *mic_u,*mic_v,*mic_w,*mic_x;
mic_u=iodata.u;
mic_v=iodata.v;
mic_w=iodata.w;
mic_x=iodata.x;
int mic_Nbase=iodata.Nbase;
int mic_tilesz=iodata.tilesz;
int mic_N=iodata.N;
double mic_freq0=iodata.freq0;
double mic_deltaf=iodata.deltaf;
double mic_data_min_uvcut=Data::min_uvcut;
int mic_data_Nt=Data::Nt;
int mic_data_max_emiter=Data::max_emiter;
int mic_data_max_iter=Data::max_iter;
int mic_data_max_lbfgs=Data::max_lbfgs;
int mic_data_lbfgs_m=Data::lbfgs_m;
int mic_data_gpu_threads=Data::gpu_threads;
int mic_data_linsolv=Data::linsolv;
int mic_data_solver_mode=Data::solver_mode;
int mic_data_randomize=Data::randomize;
double mic_data_nulow=Data::nulow;
double mic_data_nuhigh=Data::nuhigh;
#endif
/* FIXME: uvmin is not needed in calibration, because its taken care of by flags */
if (!Data::DoSim) {
@ -637,6 +708,12 @@ beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,bea
}
/****************** end calibration **************************/
/****************** begin diagnostics ************************/
#ifdef HAVE_CUDA
if (Data::DoDiag) {
calculate_diagnostics(iodata.u,iodata.v,iodata.w,p,iodata.xo,iodata.N,iodata.Nbase,iodata.tilesz,barr,carr,coh,M,Mt,Data::DoDiag,Data::Nt);
}
#endif /* HAVE_CUDA */
/****************** end diagnostics **************************/
} else {
/************ simulation only mode ***************************/
if (!solfile) {
@ -684,11 +761,33 @@ beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,bea
}
/**********************************************************/
/* also write back */
if (iodata.Nms==1) {
Data::writeData(msitr[0]->table(),iodata);
} else {
Data::writeDataList(msitr,iodata);
}
for(int cm=0; cm<iodata.Nms; cm++) {
(*msitr[cm])++;
}
if (!Data::DoSim) {
/* if residual has increased too much, or all are flagged (0 residual)
or NaN
reset solutions to original
initial values */
if (res_1==0.0 || !isfinite(res_1) || res_1>res_ratio*res_prev) {
cout<<"Resetting Solution"<<endl;
/* reset solutions so next iteration has default initial values */
memcpy(p,pinit,(size_t)iodata.N*8*Mt*sizeof(double));
/* also assume iterations have restarted from scratch */
start_iter=1;
/* also forget min residual (otherwise will try to reset it always) */
res_prev=res_1;
} else if (res_1<res_prev) { /* only store the min value */
res_prev=res_1;
}
}
end_time = time(0);
elapsed_time = ((double) (end_time-start_time)) / 60.0;
if (!Data::DoSim) {
if (solver_mode==SM_OSLM_OSRLM_RLBFGS||solver_mode==SM_RLM_RLBFGS||solver_mode==SM_RTR_OSRLM_RLBFGS || solver_mode==SM_NSD_RLBFGS) {
@ -698,6 +797,18 @@ beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,bea
} else {
cout<<"Timeslot: "<<tilex<<", Time spent="<<elapsed_time<<" minutes"<<endl;
}
#ifdef HAVE_CUDA
/* if -E uses a large value ~say 100, at each multiple of this, clear GPU memory */
if (GPUpredict>1 && tilex>0 && !(tilex%GPUpredict)) {
for (int gpuid=0; gpuid<=MAX_GPU_ID; gpuid++) {
cudaSetDevice(gpuid);
cudaDeviceReset();
cudaDeviceSetLimit(cudaLimitMallocHeapSize, Data::heapsize*1024*1024);
}
}
#endif
}
@ -712,6 +823,10 @@ beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,bea
Data::freeData(iodata,beam);
}
#ifdef USE_MIC
free(mic_pindex);
free(mic_chunks);
#endif
/**********************************************************/
exinfo_gaussian *exg;
@ -784,9 +899,6 @@ beam.p_ra0,beam.p_dec0,iodata.freq0,beam.sx,beam.sy,beam.time_utc,beam.Nelem,bea
}
/**********************************************************/
#ifdef HAVE_CUDA
cudaDeviceReset();
#endif
cout<<"Done."<<endl;
return 0;
}