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

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

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