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

Last change on this file since 3892 was 2328, checked in by aubourg, 23 years ago

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