source: Sophya/trunk/ArchTOIPipe/TestPipes/quickmap_p.cc@ 1826

Last change on this file since 1826 was 1820, checked in by aubourg, 24 years ago

GPH 424
File size: 12.7 KB
Line 
1// GPH 424.1 Planck HFI-L2 Simple Map Making
2// Parallel version
3// Eric Aubourg CEA/DAPNIA/SPP
4// $Id: quickmap_p.cc,v 1.4 2001-12-17 23:33:46 aubourg Exp $
5
6#include <stdlib.h>
7#include <stdio.h>
8#include "toi.h"
9#include "toiprocessor.h"
10#include "fitstoirdr.h"
11#include "fitstoiwtr.h"
12#include "toimanager.h"
13#include "toisegment.h"
14
15#include "sambainit.h"
16#include "toi2map.h"
17#include "fitsspherehealpix.h"
18#include "localmap.h"
19#include "fitslocalmap.h"
20
21#include <util/PlatformUtils.hpp>
22#include <sax/SAXException.hpp>
23#include <sax/SAXParseException.hpp>
24#include <parsers/DOMParser.hpp>
25#include <dom/DOM_DOMException.hpp>
26
27
28void usage(void);
29
30void usage(void) {
31 cerr << "usage: gph424_1p xml_options_file" << endl;
32 exit(-1);
33}
34
35// Prefs structure
36
37struct Gph424_1_Prefs {
38 struct {
39 struct {
40 string fname;
41 string coord1; // TOI name
42 string coord2; // TOI name
43 unsigned long typcoord;
44 } pointing;
45 struct {
46 string fname;
47 string toiname;
48 } signal;
49 } input;
50
51 struct {
52 string fname;
53 string wfname; // weight
54 unsigned long typcoord;
55 string typmap; // local, healpix...
56 struct {
57 int nlat;
58 } healpix_opt;
59 struct {
60 int npx, npy;
61 double org1,org2;
62 double anglex, angley;
63 } local_opt;
64 } output;
65
66 int nproc;
67 long snb, sne;
68
69 void dump() {
70 cout << "*** GPH424-1 settings\n";
71 cout << " input\n";
72 cout << " | pointing\n";
73 cout << " | | fname " << input.pointing.fname << "\n";
74 cout << " | | coord1 " << input.pointing.coord1 << "\n";
75 cout << " | | coord2 " << input.pointing.coord2 << "\n";
76 cout << " | | typcoord "
77 << ((input.pointing.typcoord & TypCoordGal) ? "G" : "E")
78 << ((input.pointing.typcoord & TypCoord1H) ? "H" :
79 (input.pointing.typcoord & TypCoord1D) ? "D" : "R")
80 << ((input.pointing.typcoord & TypCoord1C) ? "C" : "L")
81 << ((input.pointing.typcoord & TypCoord2H) ? "H" :
82 (input.pointing.typcoord & TypCoord2D) ? "D" : "R")
83 << ((input.pointing.typcoord & TypCoord2C) ? "C" : "L")
84 << "\n";
85 cout << " | signal\n";
86 cout << " | | fname " << input.signal.fname << "\n";
87 cout << " | | toiname " << input.signal.toiname << "\n";
88 cout << " output\n";
89 cout << " | fname " << output.fname << "\n";
90 cout << " | wfname " << output.wfname << "\n";
91 cout << " | typcoord "
92 << ((output.typcoord & TypCoordGal) ? "G" : "E")
93 << "\n";
94 cout << " | typmap " << output.typmap << "\n";
95 if (output.typmap == "healpix") {
96 cout << " | | nlat " << output.healpix_opt.nlat << "\n";
97 } else if (output.typmap == "local") {
98 cout << " | | npix " << output.local_opt.npx << ", " << output.local_opt.npy << "\n";
99 cout << " | | origin " << output.local_opt.org1 << ", " << output.local_opt.org2 << "\n";
100 cout << " | | extension " << output.local_opt.anglex << ", " << output.local_opt.angley << "\n";
101 }
102 cout << " nproc " << nproc << "\n";
103 if (sne >= 0 && snb >= 0) {
104 cout << " samples " << snb << " " << sne << "\n";
105 }
106 cout << "***" << endl;
107 }
108};
109
110DOM_Node DOMGetNode(DOM_Document& doc, string path) {
111 DOM_Element el = doc.getDocumentElement();
112
113 while (path != "") {
114 if (path[0] == '/') {
115 path = path.substr(1);
116 }
117 if (path == "") break;
118 int x = path.find('/');
119 string pathElem;
120 if (x == -1) {
121 pathElem = path;
122 path = "";
123 } else {
124 pathElem = path.substr(0, x);
125 path = path.substr(x);
126 }
127 DOM_Node n = el.getElementsByTagName(pathElem.c_str()).item(0);
128 if (n==0) return n;
129 el = (DOM_Element&) n;
130 }
131 return el;
132}
133
134string DOMGetString(DOM_Document& doc, string path, string def="") {
135 DOM_Node n = DOMGetNode(doc, path);
136 if (n==0) return def;
137 char* s = n.getFirstChild().getNodeValue().transcode();
138 char* p=s; while ((*p)==' ') p++;
139 char* q=p+strlen(p)-1;
140 while (q>p && *q==' ') {
141 *q = '\0';
142 q--;
143 }
144 string ss = p;
145 free(s);
146 return ss;
147}
148
149bool DOMHasOption(DOM_Document& doc, string path) {
150 DOM_Node n = DOMGetNode(doc, path);
151 return n != 0;
152}
153
154void parsePrefs(string prefFile, Gph424_1_Prefs& prefs) {
155 try {
156 XMLPlatformUtils::Initialize();
157 } catch (const XMLException& toCatch) {
158 cerr << "Error during XML initialization! :\n"
159 << toCatch.getMessage() << endl;
160 exit(-1);
161 }
162
163 DOMParser* parser = new DOMParser;
164
165 try {
166 parser->parse(prefFile.c_str());
167 } catch (const XMLException& toCatch) {
168 cerr << "\nError during parsing: '" << prefFile << "'\n"
169 << "Exception message is: \n"
170 << toCatch.getMessage() << "\n" << endl;
171 exit(-1);
172 } catch (const DOM_DOMException& toCatch) {
173 cerr << "\nDOM Error during parsing: '" << prefFile << "'\n"
174 << "DOMException code is: \n"
175 << toCatch.code << "\n" << endl;
176 exit(-1);
177 } catch (...) {
178 cerr << "\nUnexpected exception during parsing: '" << prefFile << "'\n";
179 exit(-1);
180 }
181
182 DOM_Document doc = parser->getDocument();
183
184 prefs.input.pointing.fname = DOMGetString(doc, "/input/pointing/fname");
185 prefs.input.pointing.coord1 = DOMGetString(doc, "/input/pointing/coord1","phi");
186 prefs.input.pointing.coord2 = DOMGetString(doc, "/input/pointing/coord2","theta");
187 bool inequ = DOMHasOption(doc, "/input/pointing/equatorial");
188 bool ingal = DOMHasOption(doc, "/input/pointing/galactic");
189 prefs.input.pointing.typcoord = ingal ? TypCoordGal : TypCoordEq;
190 bool inhd = DOMHasOption(doc, "/input/pointing/hourdeg");
191 bool indd = DOMHasOption(doc, "/input/pointing/degdeg");
192 bool inrr = DOMHasOption(doc, "/input/pointing/radian");
193 bool iscolat = DOMHasOption(doc, "/input/pointing/colat");
194 if (inhd) {
195 prefs.input.pointing.typcoord |= TypCoord1H|TypCoord2D;
196 } else if (indd) {
197 prefs.input.pointing.typcoord |= TypCoord1D|TypCoord2D;
198 } else if (inrr) {
199 prefs.input.pointing.typcoord |= TypCoord1R|TypCoord2R;
200 } else if (prefs.input.pointing.typcoord & TypCoordEq) {
201 prefs.input.pointing.typcoord = TypCoordEqStd;
202 } else {
203 prefs.input.pointing.typcoord = TypCoordGalStd;
204 }
205 if (iscolat) {
206 prefs.input.pointing.typcoord |= TypCoord2C;
207 }
208
209 prefs.input.signal.fname = DOMGetString(doc, "/input/signal/fname");
210 prefs.input.signal.toiname = DOMGetString(doc, "/input/signal/toiname");
211
212 prefs.output.fname = DOMGetString(doc, "/output/fname");
213 prefs.output.wfname = DOMGetString(doc, "/output/wfname");
214 bool outequ = DOMHasOption(doc, "/output/equatorial");
215 bool outgal = DOMHasOption(doc, "/output/galactic");
216 prefs.output.typcoord = outgal ? TypCoordGal : TypCoordEq;
217 bool outhd = DOMHasOption(doc, "/output/hourdeg");
218 bool outdd = DOMHasOption(doc, "/output/degdeg");
219 bool outrr = DOMHasOption(doc, "/output/radian");
220 if (outhd) {
221 prefs.output.typcoord |= TypCoord1H | TypCoord2D;
222 } else if (outdd) {
223 prefs.output.typcoord |= TypCoord1D | TypCoord2D;
224 } else if (outrr) {
225 prefs.output.typcoord |= TypCoord1R | TypCoord2R;
226 } else if (prefs.output.typcoord & TypCoordEq) {
227 prefs.output.typcoord = TypCoordEqStd;
228 } else {
229 prefs.output.typcoord |= TypCoordGalStd;
230 }
231
232 prefs.output.typmap = DOMHasOption(doc, "/output/local") ? "local" :
233 (DOMHasOption(doc, "/output/healpix") ? "healpix" : "healpix");
234 prefs.output.healpix_opt.nlat = atoi(DOMGetString(doc, "/output/healpix/nlat", "128").c_str());
235
236 string s = DOMGetString(doc, "/output/local/npx");
237 prefs.output.local_opt.npx = prefs.output.local_opt.npy = -1;
238 int nc = sscanf(s.c_str(), "%d %d", &prefs.output.local_opt.npx, &prefs.output.local_opt.npy);
239 if (nc == 1) prefs.output.local_opt.npy = prefs.output.local_opt.npx;
240
241 s = DOMGetString(doc, "/output/local/center");
242 prefs.output.local_opt.org1 = prefs.output.local_opt.org2 = 0;
243 sscanf(s.c_str(), "%lg %lg", &prefs.output.local_opt.org1, &prefs.output.local_opt.org2);
244
245 s = DOMGetString(doc, "/output/local/angle");
246 prefs.output.local_opt.anglex = prefs.output.local_opt.angley = 0;
247 nc = sscanf(s.c_str(), "%d %d", &prefs.output.local_opt.anglex, &prefs.output.local_opt.angley);
248 if (nc == 1) prefs.output.local_opt.angley = prefs.output.local_opt.anglex;
249
250 prefs.nproc = atoi(DOMGetString(doc, "/nproc", "1").c_str());
251
252 s = DOMGetString(doc, "/samples", "-1 -1");
253 sscanf(s.c_str(), "%ld %ld", &prefs.snb, &prefs.sne);
254}
255
256int main(int argc, char** argv) {
257
258 // Read prefs
259 Gph424_1_Prefs prefs;
260 parsePrefs(argv[1], prefs);
261 prefs.dump();
262
263 int np = prefs.nproc;
264
265 // Assemble pipeline
266
267 TOIManager* mgr = TOIManager::getManager();
268
269 // Two input files for each data chunk
270
271 vector<FITSTOIReader*> inPoint;
272 vector<FITSTOIReader*> inSig;
273 for (int i=0; i<np; i++) {
274 inPoint.push_back(new FITSTOIReader(prefs.input.pointing.fname));
275 inSig.push_back(new FITSTOIReader(prefs.input.signal.fname));
276 inPoint[i]->setImplicitSN();
277 inSig[i]->setImplicitSN();
278 }
279
280 // Prepare maps for output
281
282 vector<PixelMap<r_8>*> map;
283 vector<PixelMap<r_8>*> wmap;
284
285 for (int i=0; i<np; i++) {
286 if (prefs.output.typmap == "local") {
287 LocalMap<r_8>* lmap =
288 new LocalMap<r_8>(prefs.output.local_opt.npx, prefs.output.local_opt.npy);
289 lmap->SetOrigin(prefs.output.local_opt.org1, prefs.output.local_opt.org2);
290 lmap->SetSize(prefs.output.local_opt.anglex, prefs.output.local_opt.angley);
291 map.push_back(lmap);
292 wmap.push_back(new LocalMap<r_8>(*lmap));
293 } else if (prefs.output.typmap == "healpix") {
294 map.push_back(new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat));
295 wmap.push_back(new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat));
296 } else {
297 cout << "typmap " << prefs.output.typmap << " unknown" << endl;
298 exit(-1);
299 }
300 }
301
302 // One map projector per data chunk
303
304 vector<TOI2Map*> toi2m;
305 {for (int i=0; i<np; i++) {
306 toi2m.push_back(new TOI2Map(map[i],wmap[i]));
307 toi2m[i]->SetEquinox(2000.);
308 toi2m[i]->SetCoorIn((TypAstroCoord) prefs.input.pointing.typcoord);
309 toi2m[i]->SetCoorMap((TypAstroCoord) prefs.output.typcoord);
310 }}
311
312
313 if (prefs.snb >= 0 && prefs.sne >= 0) {
314 mgr->setRequestedSample(prefs.snb, prefs.sne);
315 }
316
317 int sn0 = inPoint[0]->getMinOut();
318 int sn1 = inPoint[0]->getMaxOut();
319
320 cout << "Range to be processed : " << sn0 << " - " << sn1 << endl;
321
322 int chunksize = (sn1-sn0)/np;
323 {for (int i=0; i<np; i++) {
324 int isn0 = sn0 + chunksize*i;
325 int isn1 = sn0 + chunksize*(i+1) - 1;
326 if (i == np-1) isn1 = sn1;
327 cout << "Processor " << i << " range " << isn0 << " - " << isn1 << endl;
328 inPoint[i]->setMinSn(isn0);
329 inSig[i]->setMinSn(isn0);
330 inPoint[i]->setMaxSn(isn1);
331 inSig[i]->setMaxSn(isn1);
332 }}
333
334 // Connect pipes
335
336 {for (int i=0; i<np; i++) {
337 TOISegmented * toicoord1in = new TOISegmented(prefs.input.pointing.coord1);
338 inPoint[i]->addOutput(prefs.input.pointing.coord1,toicoord1in);
339 toi2m[i]->addInput("Coord1In",toicoord1in);
340
341 TOISegmented * toicoord2in = new TOISegmented(prefs.input.pointing.coord2);
342 inPoint[i]->addOutput(prefs.input.pointing.coord2,toicoord2in);
343 toi2m[i]->addInput("Coord2In",toicoord2in);
344
345 TOISegmented * toibolin = new TOISegmented("toi_bolo_in");
346 inSig[i]->addOutput(prefs.input.signal.toiname,toibolin);
347 toi2m[i]->addInput("BoloIn",toibolin);
348 }}
349
350 // Start processing
351
352 {for (int i=0; i<np; i++) {
353 inPoint[i]->start();
354 inSig[i]->start();
355 toi2m[i]->start();
356 }}
357 mgr->joinAll();
358
359 // Merge maps
360
361 PixelMap<r_8>* smap = map[0];
362 PixelMap<r_8>* swmap = wmap[0];
363
364 {for (int i=1; i<np; i++) {
365 // There should be a += operator in all maps, not only in local maps
366 // *smap += *map[i];
367 // *swmap += *wmap[i];
368 int n = smap->NbPixels();
369 for (int j=0; j<n; j++) {
370 (*smap)(j) += (*map[i])(j);
371 (*swmap)(j) += (*wmap[i])(j);
372 }
373 }}
374
375 // Save maps
376
377 {
378 FitsOutFile sfits(prefs.output.fname,FitsFile::clear);
379 cout<<"Creating map fits file "<<prefs.output.fname<<endl;
380 // Would need some generic mechanism here $CHECK$ $$TBD$$
381 if ((string)smap->TypeOfMap() == "RING") {
382 sfits << * (SphereHEALPix<r_8>*)smap;
383 } else if ((string)smap->TypeOfMap() == "LOCAL") {
384 sfits << * (LocalMap<r_8>*)smap;
385 } else {
386 cout << "error, unknown type " << smap->TypeOfMap() << endl;
387 }
388 }
389
390 {
391 FitsOutFile swfits(prefs.output.wfname,FitsFile::clear);
392 cout<<"Creating map weight fits file "<<prefs.output.wfname<<endl;
393 // swfits << *wmap;
394 if ((string)swmap->TypeOfMap() == "RING") {
395 swfits << * (SphereHEALPix<r_8>*)swmap;
396 } else if ((string)swmap->TypeOfMap() == "LOCAL") {
397 swfits << * (LocalMap<r_8>*)swmap;
398 } else {
399 cout << "error, unknown type " << swmap->TypeOfMap() << endl;
400 }
401 }
402
403 return 0;
404}
Note: See TracBrowser for help on using the repository browser.