Healpix C++  3.50
hpxtest.cc
1 /*
2  * This file is part of Healpix_cxx.
3  *
4  * Healpix_cxx is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * Healpix_cxx is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with Healpix_cxx; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  *
18  * For more information about HEALPix, see http://healpix.sourceforge.net
19  */
20 
21 /*
22  * Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
23  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
24  * (DLR).
25  */
26 
27 /*
28  * Copyright (C) 2004-2017 Max-Planck-Society
29  * Author: Martin Reinecke
30  */
31 
32 /*
33 
34 Candidates for testing the validity of the Healpix routines:
35 
36 - done: ang2pix(pix2ang(i)) = i for all Healpix_Bases
37 - done: pix2ang(ang2pix(ptg)) dot ptg > 1-delta for all Healpix_Bases
38 - done: ring2nest(nest2ring(i)) = i for all hierarchical Healpix_Bases
39 - done: downgrade(upgrade(map)) = map for all maps
40 - done: map and downgraded map should have same average
41 - done: alm2map(map2alm(map)) approx map (same for pol)
42 - partly done: neighbor tests
43 - powspec -> alm -> powspec (should produce similar powspecs, also for pol)
44 - done: two swap_schemes() should produce original map
45 - done: query_disc tests (dot products etc.)
46 - a_lms: test Set(), Scale(), Add(), alm(l,m) = alm.mstart(m)[l], etc.
47 
48 */
49 
50 #include <iostream>
51 #include "healpix_base.h"
52 #include "healpix_map.h"
53 #include "arr.h"
54 #include "planck_rng.h"
55 #include "lsconstants.h"
56 #include "alm.h"
57 #include "alm_healpix_tools.h"
58 #include "alm_powspec_tools.h"
59 #include "geom_utils.h"
60 #include "walltimer.h"
61 #include "announce.h"
62 #include "compress_utils.h"
63 #include "moc.h"
64 #include "moc_fitsio.h"
65 #include "crangeset.h"
66 #include "weight_utils.h"
67 #include "powspec.h"
68 
69 using namespace std;
70 
71 #define UNITTESTS
72 
73 #ifdef UNITTESTS
74 int errcount=0;
75 #define FAIL(a) {a; if (++errcount>100) planck_fail("unit test errors");}
76 #else
77 #define FAIL(a) a;
78 #endif
79 
80 const int nsamples = 1000000;
81 
82 planck_rng rng;
83 
84 namespace {
85 
86 void random_dir (pointing &ptg)
87  {
88  ptg.theta = acos(rng.rand_uni()*2-1);
89  ptg.phi = rng.rand_uni()*twopi;
90  }
91 void random_zphi (double &z, double &phi)
92  {
93  z = rng.rand_uni()*2-1;
94  phi = rng.rand_uni()*twopi;
95  }
96 
97 template<typename I> string bname()
98  { return string("(basetype: ")+type2typename<I>()+")"; }
99 
100 
101 template<typename I> rangeset<I> randomRangeSet(int num, I start, int dist)
102  {
103  rangeset<I> rs;
104  I curval=start;
105  for (int i=0; i<num; ++i)
106  {
107  I v1=curval+1+(rng.int_rand_uni()%dist);
108  I v2=v1+1+(rng.int_rand_uni()%dist);
109  rs.append(v1,v2);
110  curval=v2;
111  }
112  return rs;
113  }
114 template<typename I> Moc<I> randomMoc(int num, I start, int dist)
115  {
116  Moc<I> moc;
117  I curval=start+(I(1)<<(2*Moc<I>::maxorder));
118  for (int i=0; i<num; ++i)
119  {
120  I v1=curval+1+(rng.int_rand_uni()%dist);
121  I v2=v1+1+(rng.int_rand_uni()%dist);
122  moc.addPixelRange(Moc<I>::maxorder,v1,v2);
123  curval=v2;
124  }
125  return moc;
126  }
127 
128 template<typename I> rangeset<I> makeRS(const string &input)
129  {
130  istringstream is(input);
131  rangeset<I> res;
132  I a,b;
133  while (is)
134  {
135  is>>a>>b;
136  planck_assert (is||is.eof(),"aargh");
137  if (is) res.append(a,b);
138  }
139  return res;
140  }
141 template<typename I> void rsOps(const rangeset<I> &a, const rangeset<I> &b)
142  {
143  planck_assert(!b.overlaps(a.op_andnot(b)),"error");
144  planck_assert(a.op_or(b).nval()>=a.nval(),"error");
145  planck_assert(a.op_or(b).nval()>=b.nval(),"error");
146  planck_assert(a.op_and(b).nval()<=a.nval(),"error");
147  planck_assert(a.op_and(b).nval()<=b.nval(),"error");
148  planck_assert(a.op_or(b).contains(a),"error");
149  planck_assert(a.op_or(b).contains(b),"error");
150  planck_assert(!a.op_andnot(b).overlaps(b),"error");
151  planck_assert(!b.op_andnot(a).overlaps(a),"error");
152  planck_assert(a.op_or(b).nval()==a.nval()+b.nval()-a.op_and(b).nval(),
153  "error");
154  planck_assert(a.op_or(b).op_andnot(a.op_and(b)).nval()==
155  a.op_or(b).nval()-a.op_and(b).nval(),"error");
156  planck_assert(a.op_xor(b)==a.op_andnot(b).op_or(b.op_andnot(a)),"error");
157  }
158 template<typename I> void check_rangeset()
159  {
160  cout << "testing rangeset " << bname<I>() << endl;
161  {
162  rangeset<I> b;
163  b.append(1,11);
164  planck_assert(b==makeRS<I>("1 11"),"mismatch");
165  b.append(10,15);
166  planck_assert(b==makeRS<I>("1 15"),"mismatch");
167  b.append(1,15);
168  planck_assert(b==makeRS<I>("1 15"),"mismatch");
169  b.append(7,15);
170  planck_assert(b==makeRS<I>("1 15"),"mismatch");
171  b.append(30,41);
172 #if 0
173  planck_assert(b==makeRS<I>("1 15 30 41"),"mismatch");
174  try
175  {
176  b.append(29,31);
177  planck_fail("Should have raised an IllegalArgumentException");
178  }
179  catch (PlanckError E) {}
180 #endif
181  }
182  {
183  rangeset<I> b=makeRS<I>("1 11 30 41");
184  planck_assert(!b.contains(0),"error");
185  planck_assert(b.contains(1),"error");
186  planck_assert(b.contains(5),"error");
187  planck_assert(b.contains(10),"error");
188  planck_assert(!b.contains(11),"error");
189  planck_assert(!b.contains(29),"error");
190  planck_assert(b.contains(30),"error");
191  planck_assert(b.contains(35),"error");
192  planck_assert(b.contains(40),"error");
193  planck_assert(!b.contains(41),"errror");
194  }
195  {
196  rangeset<I> b;
197  b.add(5, 11);
198  planck_assert(b==makeRS<I>("5 11"),"error");
199  b.add(1, 7);
200  planck_assert(b==makeRS<I>("1 11"),"error");
201  b.add(1, 11);
202  planck_assert(b==makeRS<I>("1 11"),"error");
203  b.add(30, 41);
204  planck_assert(b==makeRS<I>("1 11 30 41"),"error");
205  b.add(1, 11);
206  planck_assert(b==makeRS<I>("1 11 30 41"),"error");
207  b.add(-1,0);
208  planck_assert(b==makeRS<I>("-1 0 1 11 30 41"),"error");
209  b.add(-2,-1);
210  planck_assert(b==makeRS<I>("-2 0 1 11 30 41"),"error");
211  b.add(-2,-1);
212  planck_assert(b==makeRS<I>("-2 0 1 11 30 41"),"error");
213  b.add(2, 11);
214  planck_assert(b==makeRS<I>("-2 0 1 11 30 41"),"error");
215  b.add(1, 10);
216  planck_assert(b==makeRS<I>("-2 0 1 11 30 41"),"error");
217  b.add(15, 21);
218  planck_assert(b==makeRS<I>("-2 0 1 11 15 21 30 41"),"error");
219  }
220  {
221  rangeset<I> b = makeRS<I>("0 11 20 31");
222  b.remove(5,25);
223  planck_assert(b==makeRS<I>("0 5 25 31"),"error");
224  b.remove(31,32);
225  planck_assert(b==makeRS<I>("0 5 25 31"),"error");
226  b.remove(35,38);
227  planck_assert(b==makeRS<I>("0 5 25 31"),"error");
228  b.remove(-90,-80);
229  planck_assert(b==makeRS<I>("0 5 25 31"),"error");
230  b.remove(27,29);
231  planck_assert(b==makeRS<I>("0 5 25 27 29 31"),"error");
232  b.remove(25,26);
233  planck_assert(b==makeRS<I>("0 5 26 27 29 31"),"error");
234  b.remove(4,6);
235  planck_assert(b==makeRS<I>("0 4 26 27 29 31"),"error");
236  b.remove(-20,40);
237  planck_assert(b==makeRS<I>(""),"error");
238  b.remove(-20,40);
239  planck_assert(b==makeRS<I>(""),"error");
240  }
241  {
242  rangeset<I> b = makeRS<I>("0 11 20 31");
243  b.intersect(2,29);
244  planck_assert(b==makeRS<I>("2 11 20 29"),"error");
245  b.intersect(-8,50);
246  planck_assert(b==makeRS<I>("2 11 20 29"),"error");
247  b.intersect(2,50);
248  planck_assert(b==makeRS<I>("2 11 20 29"),"error");
249  b.intersect(2,29);
250  planck_assert(b==makeRS<I>("2 11 20 29"),"error");
251  b.intersect(-18,29);
252  planck_assert(b==makeRS<I>("2 11 20 29"),"error");
253  b.intersect(3,11);
254  planck_assert(b==makeRS<I>("3 11"),"error");
255  b = makeRS<I>("0 11 20 31");
256  b.intersect(3,15);
257  planck_assert(b==makeRS<I>("3 11"),"error");
258  b = makeRS<I>("0 11 20 31");
259  b.intersect(17,30);
260  planck_assert(b==makeRS<I>("20 30"),"error");
261  b = makeRS<I>("0 11 20 31");
262  b.intersect(11,20);
263  planck_assert(b==makeRS<I>(""),"error");
264  b = makeRS<I>("0 11 20 31");
265  b.intersect(-8,-7);
266  planck_assert(b==makeRS<I>(""),"error");
267  b = makeRS<I>("0 11 20 31");
268  b.intersect(31,35);
269  planck_assert(b==makeRS<I>(""),"error");
270  }
271  {
272  planck_assert(makeRS<I>("1 11 20 31 40 56")==
273  makeRS<I>("20 31 40 51").op_or
274  (makeRS<I>("1 11 45 56")),"error");
275  planck_assert(makeRS<I>("1 11 45 56")==
276  makeRS<I>("").op_or
277  (makeRS<I>("1 11 45 56")),"error");
278  planck_assert(makeRS<I>("1 11 45 56")==
279  makeRS<I>("1 11 45 56").op_or
280  (makeRS<I>("")),"error");
281  }
282  {
283  planck_assert(makeRS<I>("22 24 45 51")==
284  makeRS<I>("20 31 40 51").op_and
285  (makeRS<I>("1 11 22 24 45 56")),"error");
286  planck_assert(makeRS<I>("20 31 40 51 90 101 110 121 200 201")==
287  makeRS<I>("10 101 110 121 200 221").op_and
288  (makeRS<I>("20 31 40 51 90 201")),"error");
289  planck_assert(makeRS<I>("")==
290  makeRS<I>("20 31 40 51").op_and
291  (makeRS<I>("")),"error");
292  planck_assert(makeRS<I>("")==
293  makeRS<I>("").op_and
294  (makeRS<I>("20 31 40 51")),"error");
295  }
296  {
297  planck_assert(makeRS<I>("20 31 40 45")==
298  makeRS<I>("20 31 40 51").op_andnot
299  (makeRS<I>("1 11 45 56")),"error");
300  planck_assert(makeRS<I>("")==
301  makeRS<I>("").op_andnot
302  (makeRS<I>("1 11 45 56")),"error");
303  planck_assert(makeRS<I>("1 11 45 56")==
304  makeRS<I>("1 11 45 56").op_andnot
305  (makeRS<I>("")),"error");
306  }
307  {
308  rangeset<I> b = makeRS<I>("20 31 40 51");
309 
310  planck_assert(!b.contains(0,11),"error");
311  planck_assert(!b.contains(10,21),"error");
312  planck_assert(!b.contains(19,20),"error");
313  planck_assert(b.contains(20,21),"error");
314  planck_assert(b.contains(21,22),"error");
315  planck_assert(b.contains(20,31),"error");
316  planck_assert(!b.contains(25,36),"error");
317  planck_assert(b.contains(30,31),"error");
318  planck_assert(!b.contains(31,32),"error");
319  planck_assert(!b.contains(35,38),"error");
320  planck_assert(!b.contains(35,46),"error");
321  planck_assert(b.contains(40,41),"error");
322  planck_assert(!b.contains(45,56),"error");
323  planck_assert(!b.contains(60,71),"error");
324  }
325  {
326  rangeset<I> b = makeRS<I>("20 31 40 51");
327 
328  planck_assert(b.contains(makeRS<I>("20 31 40 51")),"error");
329  planck_assert(b.contains(makeRS<I>("20 21")),"error");
330  planck_assert(b.contains(makeRS<I>("50 51")),"error");
331  planck_assert(!b.contains(makeRS<I>("19 31 40 51")),"error");
332  planck_assert(!b.contains(makeRS<I>("20 31 40 52")),"error");
333  planck_assert(!b.contains(makeRS<I>("20 51")),"error");
334  planck_assert(!b.contains(makeRS<I>("0 1")),"error");
335  planck_assert(!b.contains(makeRS<I>("0 20 31 40 51 100")),"error");
336  }
337  {
338  rangeset<I> b = makeRS<I>("20 31 40 51");
339 
340  planck_assert(!b.overlaps(0,11),"error");
341  planck_assert(b.overlaps(10,21),"error");
342  planck_assert(!b.overlaps(19,20),"error");
343  planck_assert(b.overlaps(20,21),"error");
344  planck_assert(b.overlaps(21,22),"error");
345  planck_assert(b.overlaps(20,31),"error");
346  planck_assert(b.overlaps(25,36),"error");
347  planck_assert(b.overlaps(30,37),"error");
348  planck_assert(!b.overlaps(31,32),"error");
349  planck_assert(!b.overlaps(35,38),"error");
350  planck_assert(b.overlaps(35,46),"error");
351  planck_assert(b.overlaps(40,41),"error");
352  planck_assert(b.overlaps(45,56),"error");
353  planck_assert(!b.overlaps(60,71),"error");
354  }
355  {
356  rangeset<I> b = makeRS<I>("20 31 40 51");
357 
358  planck_assert(b.overlaps(makeRS<I>("20 31 40 51")),"error");
359  planck_assert(b.overlaps(makeRS<I>("20 21")),"error");
360  planck_assert(b.overlaps(makeRS<I>("50 51")),"error");
361  planck_assert(b.overlaps(makeRS<I>("19 31 40 51")),"error");
362  planck_assert(b.overlaps(makeRS<I>("20 31 40 52")),"error");
363  planck_assert(b.overlaps(makeRS<I>("20 51")),"error");
364  planck_assert(!b.overlaps(makeRS<I>("0 1")),"error");
365  planck_assert(!b.overlaps(makeRS<I>("0 20 31 40 51 100")),"error");
366  }
367  {
368  const int niter = 1000;
369  for (int iter=0; iter<niter; ++iter)
370  {
371  rangeset<I> a = randomRangeSet<I>(1000, 0, 100);
372  rangeset<I> b = randomRangeSet<I>(1000, 0, 100);
373  rangeset<I> c = randomRangeSet<I>(10, 10000, 100);
374  rsOps(a,b);
375  rsOps(a,c);
376  rsOps(c,a);
377  }
378  }
379  }
380 template<typename I> crangeset<I> randomCRangeSet(int num, I start, int dist)
381  {
382  crangeset<I> rs;
383  I curval=start;
384  for (int i=0; i<num; ++i)
385  {
386  I v1=curval+1+(rng.int_rand_uni()%dist);
387  I v2=v1+1+(rng.int_rand_uni()%dist);
388  rs.append(v1,v2);
389  curval=v2;
390  }
391  return rs;
392  }
393 template<typename I> crangeset<I> makeCRS(const string &input)
394  {
395  istringstream is(input);
396  crangeset<I> res;
397  I a,b;
398  while (is)
399  {
400  is>>a>>b;
401  planck_assert (is||is.eof(),"aargh");
402  if (is) res.append(a,b);
403  }
404  return res;
405  }
406 template<typename I> void crsOps(const crangeset<I> &a, const crangeset<I> &b)
407  {
408  planck_assert(!b.overlaps(a.op_andnot(b)),"error");
409  planck_assert(a.op_or(b).nval()>=a.nval(),"error");
410  planck_assert(a.op_or(b).nval()>=b.nval(),"error");
411  planck_assert(a.op_and(b).nval()<=a.nval(),"error");
412  planck_assert(a.op_and(b).nval()<=b.nval(),"error");
413  planck_assert(a.op_or(b).contains(a),"error");
414  planck_assert(a.op_or(b).contains(b),"error");
415  planck_assert(!a.op_andnot(b).overlaps(b),"error");
416  planck_assert(!b.op_andnot(a).overlaps(a),"error");
417  planck_assert(a.op_or(b).nval()==a.nval()+b.nval()-a.op_and(b).nval(),
418  "error");
419  planck_assert(a.op_or(b).op_andnot(a.op_and(b)).nval()==
420  a.op_or(b).nval()-a.op_and(b).nval(),"error");
421  planck_assert(a.op_xor(b)==a.op_andnot(b).op_or(b.op_andnot(a)),"error");
422  }
423 template<typename I> void check_crangeset()
424  {
425  cout << "testing crangeset " << bname<I>() << endl;
426  {
427  crangeset<I> b;
428  b.append(1,11);
429  planck_assert(b==makeCRS<I>("1 11"),"mismatch");
430  b.append(10,15);
431  planck_assert(b==makeCRS<I>("1 15"),"mismatch");
432  b.append(1,15);
433  planck_assert(b==makeCRS<I>("1 15"),"mismatch");
434  b.append(7,15);
435  planck_assert(b==makeCRS<I>("1 15"),"mismatch");
436  b.append(30,41);
437 #if 0
438  planck_assert(b==makeRS<I>("1 15 30 41"),"mismatch");
439  try
440  {
441  b.append(29,31);
442  planck_fail("Should have raised an IllegalArgumentException");
443  }
444  catch (PlanckError E) {}
445 #endif
446  }
447  {
448  crangeset<I> b=makeCRS<I>("1 11 30 41");
449  planck_assert(!b.contains(0),"error");
450  planck_assert(b.contains(1),"error");
451  planck_assert(b.contains(5),"error");
452  planck_assert(b.contains(10),"error");
453  planck_assert(!b.contains(11),"error");
454  planck_assert(!b.contains(29),"error");
455  planck_assert(b.contains(30),"error");
456  planck_assert(b.contains(35),"error");
457  planck_assert(b.contains(40),"error");
458  planck_assert(!b.contains(41),"errror");
459  }
460  {
461  planck_assert(makeCRS<I>("1 11 20 31 40 56")==
462  makeCRS<I>("20 31 40 51").op_or
463  (makeCRS<I>("1 11 45 56")),"error");
464  planck_assert(makeCRS<I>("1 11 45 56")==
465  makeCRS<I>("").op_or
466  (makeCRS<I>("1 11 45 56")),"error");
467  planck_assert(makeCRS<I>("1 11 45 56")==
468  makeCRS<I>("1 11 45 56").op_or
469  (makeCRS<I>("")),"error");
470  }
471  {
472  planck_assert(makeCRS<I>("22 24 45 51")==
473  makeCRS<I>("20 31 40 51").op_and
474  (makeCRS<I>("1 11 22 24 45 56")),"error");
475  planck_assert(makeCRS<I>("20 31 40 51 90 101 110 121 200 201")==
476  makeCRS<I>("10 101 110 121 200 221").op_and
477  (makeCRS<I>("20 31 40 51 90 201")),"error");
478  planck_assert(makeCRS<I>("")==
479  makeCRS<I>("20 31 40 51").op_and
480  (makeCRS<I>("")),"error");
481  planck_assert(makeCRS<I>("")==
482  makeCRS<I>("").op_and
483  (makeCRS<I>("20 31 40 51")),"error");
484  }
485  {
486  planck_assert(makeCRS<I>("20 31 40 45")==
487  makeCRS<I>("20 31 40 51").op_andnot
488  (makeCRS<I>("1 11 45 56")),"error");
489  planck_assert(makeCRS<I>("")==
490  makeCRS<I>("").op_andnot
491  (makeCRS<I>("1 11 45 56")),"error");
492  planck_assert(makeCRS<I>("1 11 45 56")==
493  makeCRS<I>("1 11 45 56").op_andnot
494  (makeCRS<I>("")),"error");
495  }
496  {
497  crangeset<I> b = makeCRS<I>("20 31 40 51");
498 
499  planck_assert(!b.contains(0,11),"error");
500  planck_assert(!b.contains(10,21),"error");
501  planck_assert(!b.contains(19,20),"error");
502  planck_assert(b.contains(20,21),"error");
503  planck_assert(b.contains(21,22),"error");
504  planck_assert(b.contains(20,31),"error");
505  planck_assert(!b.contains(25,36),"error");
506  planck_assert(b.contains(30,31),"error");
507  planck_assert(!b.contains(31,32),"error");
508  planck_assert(!b.contains(35,38),"error");
509  planck_assert(!b.contains(35,46),"error");
510  planck_assert(b.contains(40,41),"error");
511  planck_assert(!b.contains(45,56),"error");
512  planck_assert(!b.contains(60,71),"error");
513  }
514  {
515  crangeset<I> b = makeCRS<I>("20 31 40 51");
516 
517  planck_assert(b.contains(makeCRS<I>("20 31 40 51")),"error");
518  planck_assert(b.contains(makeCRS<I>("20 21")),"error");
519  planck_assert(b.contains(makeCRS<I>("50 51")),"error");
520  planck_assert(!b.contains(makeCRS<I>("19 31 40 51")),"error");
521  planck_assert(!b.contains(makeCRS<I>("20 31 40 52")),"error");
522  planck_assert(!b.contains(makeCRS<I>("20 51")),"error");
523  planck_assert(!b.contains(makeCRS<I>("0 1")),"error");
524  planck_assert(!b.contains(makeCRS<I>("0 20 31 40 51 100")),"error");
525  }
526  {
527  crangeset<I> b = makeCRS<I>("20 31 40 51");
528 
529  planck_assert(!b.overlaps(0,11),"error");
530  planck_assert(b.overlaps(10,21),"error");
531  planck_assert(!b.overlaps(19,20),"error");
532  planck_assert(b.overlaps(20,21),"error");
533  planck_assert(b.overlaps(21,22),"error");
534  planck_assert(b.overlaps(20,31),"error");
535  planck_assert(b.overlaps(25,36),"error");
536  planck_assert(b.overlaps(30,37),"error");
537  planck_assert(!b.overlaps(31,32),"error");
538  planck_assert(!b.overlaps(35,38),"error");
539  planck_assert(b.overlaps(35,46),"error");
540  planck_assert(b.overlaps(40,41),"error");
541  planck_assert(b.overlaps(45,56),"error");
542  planck_assert(!b.overlaps(60,71),"error");
543  }
544  {
545  crangeset<I> b = makeCRS<I>("20 31 40 51");
546 
547  planck_assert(b.overlaps(makeCRS<I>("20 31 40 51")),"error");
548  planck_assert(b.overlaps(makeCRS<I>("20 21")),"error");
549  planck_assert(b.overlaps(makeCRS<I>("50 51")),"error");
550  planck_assert(b.overlaps(makeCRS<I>("19 31 40 51")),"error");
551  planck_assert(b.overlaps(makeCRS<I>("20 31 40 52")),"error");
552  planck_assert(b.overlaps(makeCRS<I>("20 51")),"error");
553  planck_assert(!b.overlaps(makeCRS<I>("0 1")),"error");
554  planck_assert(!b.overlaps(makeCRS<I>("0 20 31 40 51 100")),"error");
555  }
556  {
557  const int niter = 1000;
558  for (int iter=0; iter<niter; ++iter)
559  {
560  crangeset<I> a = randomCRangeSet<I>(1000, 0, 100);
561  crangeset<I> b = randomCRangeSet<I>(1000, 0, 100);
562  crangeset<I> c = randomCRangeSet<I>(10, 10000, 100);
563  crsOps(a,b);
564  crsOps(a,c);
565  crsOps(c,a);
566  }
567  }
568  }
569 template<typename I> void check_Moc()
570  {
571  cout << "testing MOC " << bname<I>() << endl;
572  Moc<I> moc;
573  moc.addPixelRange(0,4,5);
574  moc.addPixelRange(0,6,7);
575  moc.addPixelRange(2,4,17);
576  moc.addPixelRange(10,3000000,3000001);
577 
578  planck_assert(moc==moc.complement().complement(),"error");
579  planck_assert(moc==Moc<I>::fromUniq(moc.toUniq()),"error");
580  planck_assert(moc.maxOrder()==10,"error");
581  Moc<I> xtmp = moc.degradedToOrder(8,false);
582  planck_assert(moc.contains(xtmp),"error");
583  planck_assert(!xtmp.contains(moc),"error");
584  planck_assert(xtmp.overlaps(moc),"error");
585  xtmp=moc.degradedToOrder(8,true);
586  planck_assert(!moc.contains(xtmp),"error");
587  planck_assert(xtmp.contains(moc),"error");
588  planck_assert(xtmp.overlaps(moc),"error");
589  planck_assert(moc==Moc<I>::fromCompressed(moc.toCompressed()),"error");
590 #if 0
591  assertEquals("inconsistency",moc,MocUtil.mocFromString(" 0/4, 6 2/ \t 4 -16 10/3000000 \t\n "));
592  assertEquals("inconsistency",moc,MocUtil.mocFromString("0/6 2/ 5 2/4 2/6- 16 0/4 10/3000000"));
593  assertEquals("inconsistency",moc,MocUtil.mocFromString
594  ("{\"0\":[6] , \"2\": [5 ], \"2\":[ 4,6,7,8,9,10,11,12,13,14,15,16], \"0\":[4], \"10\":[3000000]}"));
595  assertEquals("inconsistency",moc,MocUtil.mocFromString(MocUtil.mocToStringASCII(moc)));
596  assertEquals("inconsistency",moc,MocUtil.mocFromString(MocUtil.mocToStringJSON(moc)));
597  ByteArrayOutputStream out= new ByteArrayOutputStream();
598  MocUtil.mocToFits(moc,out);
599  ByteArrayInputStream inp = new ByteArrayInputStream(out.toByteArray());
600  assertEquals("inconsistency",moc,MocUtil.mocFromFits(inp));
601 #endif
602  {
603  tsize niter = 100;
604  Moc<I> full; full.addPixelRange(0,0,12);
605  Moc<I> empty;
606  for (tsize iter=0; iter<niter; ++iter)
607  {
608  Moc<I> a = randomMoc<I>(1000, 0, 100);
609  planck_assert(a.complement().complement()==a,"error");
610  planck_assert(!a.overlaps(a.complement()),"error");
611  planck_assert(a.op_or(a.complement())==full,"error");
612  planck_assert(a.op_and(a.complement())==empty,"error");
613 #if 0
614  write_Moc_to_fits("!healpixtestmoctmp",a);
615  planck_assert(a==read_Moc_from_fits<I>("healpixtestmoctmp"),"FITS problem");
616 #endif
617  }
618  }
619  }
620 template<typename I> void check_compress()
621  {
622  cout << "testing interpolation coding " << bname<I>() << endl;
623  planck_assert(trailingZeros(4)==2,"error");
624  planck_assert(trailingZeros(5)==0,"error");
625  planck_assert(trailingZeros(int64(1)<<48)==48,"error");
626  for (tsize x=0; x<100; ++x)
627  {
628  rangeset<I> a = randomRangeSet<I>(1000, 0, 100);
629  rangeset<I> b;
630  for (tsize i=0; i<a.nranges(); ++i)
631  b.append(a.ivbegin(i)<<6,a.ivend(i)<<6);
632  const vector<I> &v=b.data();
633  obitstream obs;
634  interpol_encode(v.begin(),v.end(),obs);
635  vector<uint8> comp=obs.state();
636  vector<I> v2;
637  ibitstream ibs(comp);
638  interpol_decode(v2,ibs);
639  planck_assert(v==v2,"data mismatch");
640  }
641  }
642 
643 template<typename I> void check_ringnestring()
644  {
645  cout << "testing ring2nest(nest2ring(m))==m " << bname<I>() << endl;
646  for (int order=0; order<=T_Healpix_Base<I>::order_max; ++order)
647  {
648  T_Healpix_Base<I> base (order,RING);
649  for (int m=0; m<nsamples; ++m)
650  {
651  I pix = I(rng.rand_uni()*base.Npix());
652  if (base.ring2nest(base.nest2ring(pix))!=pix)
653  FAIL(cout<<" PROBLEM: order = "<<order<<", pixel = "<<pix<<endl)
654  }
655  }
656  }
657 
658 template<typename I> void check_nestpeanonest()
659  {
660  cout << "testing peano2nest(nest2peano(m))==m " << bname<I>() << endl;
661  for (int order=0; order<=T_Healpix_Base<I>::order_max; ++order)
662  {
663  T_Healpix_Base<I> base (order,NEST);
664  for (int m=0; m<nsamples; ++m)
665  {
666  I pix = I(rng.rand_uni()*base.Npix());
667  if (base.peano2nest(base.nest2peano(pix))!=pix)
668  FAIL(cout<<" PROBLEM: order = "<<order<<", pixel = "<<pix<<endl)
669  }
670  }
671  }
672 
673 template<typename I> void check_pixzphipix()
674  {
675  cout << "testing zphi2pix(pix2zphi(m))==m " << bname<I>() << endl;
677  for (int order=0; order<=omax; ++order)
678  {
679  T_Healpix_Base<I> base1 (order,RING), base2 (order,NEST);
680  for (int m=0; m<nsamples; ++m)
681  {
682  double z,phi;
683  I pix = I(rng.rand_uni()*base1.Npix());
684  base1.pix2zphi(pix,z,phi);
685  if (base1.zphi2pix(z,phi)!=pix)
686  FAIL(cout<<" PROBLEM: order = "<<order<<", pixel = "<<pix<<endl)
687  base2.pix2zphi(pix,z,phi);
688  if (base2.zphi2pix(z,phi)!=pix)
689  FAIL(cout<<" PROBLEM: order = "<<order<<", pixel = "<<pix<<endl)
690  }
691  }
692  for (I nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
693  {
694  T_Healpix_Base<I> base (nside,RING,SET_NSIDE);
695  for (int m=0; m<nsamples; ++m)
696  {
697  double z,phi;
698  I pix = I(rng.rand_uni()*base.Npix());
699  base.pix2zphi(pix,z,phi);
700  if (base.zphi2pix(z,phi)!=pix)
701  FAIL(cout<<" PROBLEM: nside = "<<nside<<", pixel = "<<pix<<endl)
702  }
703  }
704  }
705 
706 template<typename I> void check_zphipixzphi()
707  {
708  cout << "testing pix2zphi(zphi2pix(ptg)) approx zphi " << bname<I>() << endl;
710  for (int order=0; order<=omax; ++order)
711  {
712  T_Healpix_Base<I> base1 (order,NEST), base2 (order,RING);
713  double mincos = min (cos(base1.max_pixrad()),0.999999999999999);
714  for (int m=0; m<nsamples; ++m)
715  {
716  double z,phi,z2,phi2;
717  random_zphi (z,phi);
718  base1.pix2zphi(base1.zphi2pix(z,phi),z2,phi2);
719  if (cosdist_zphi(z,phi,z2,phi2)<mincos)
720  FAIL(cout << " PROBLEM: order = " << order
721  << ", zphi = " << z << ", " << phi << endl)
722  base2.pix2zphi(base2.zphi2pix(z,phi),z2,phi2);
723  if (cosdist_zphi(z,phi,z2,phi2)<mincos)
724  FAIL(cout << " PROBLEM: order = " << order
725  << ", zphi = " << z << ", " << phi << endl)
726  }
727  }
728  for (int nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
729  {
730  T_Healpix_Base<I> base (nside,RING,SET_NSIDE);
731  double mincos = min (cos(base.max_pixrad()),0.999999999999999);
732  for (int m=0; m<nsamples; ++m)
733  {
734  double z,phi,z2,phi2;
735  random_zphi (z,phi);
736  base.pix2zphi(base.zphi2pix(z,phi),z2,phi2);
737  if (cosdist_zphi(z,phi,z2,phi2)<mincos)
738  FAIL(cout << " PROBLEM: nside = " << nside
739  << ", zphi = " << z << ", " << phi << endl)
740  }
741  }
742  }
743 
744 template<typename I> void check_pixangpix()
745  {
746  cout << "testing ang2pix(pix2ang(m))==m " << bname<I>() << endl;
748  for (int order=0; order<=omax; ++order)
749  {
750  T_Healpix_Base<I> base1 (order,RING), base2 (order,NEST);
751  for (int m=0; m<nsamples; ++m)
752  {
753  I pix = I(rng.rand_uni()*base1.Npix());
754  if (base1.ang2pix(base1.pix2ang(pix))!=pix)
755  FAIL(cout<<" PROBLEM: order = "<<order<<", pixel = "<<pix<<endl)
756  if (base2.ang2pix(base2.pix2ang(pix))!=pix)
757  FAIL(cout<<" PROBLEM: order = "<<order<<", pixel = "<<pix<<endl)
758  }
759  }
760  for (I nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
761  {
762  T_Healpix_Base<I> base (nside,RING,SET_NSIDE);
763  for (int m=0; m<nsamples; ++m)
764  {
765  I pix = I(rng.rand_uni()*base.Npix());
766  if (base.ang2pix(base.pix2ang(pix))!=pix)
767  FAIL(cout<<" PROBLEM: nside = "<<nside<<", pixel = "<<pix<<endl)
768  }
769  }
770  }
771 
772 template<typename I> void check_neighbors()
773  {
774  cout << "testing neighbor function " << bname<I>() << endl;
776  for (int order=0; order<=omax; ++order)
777  {
778  T_Healpix_Base<I> base (order,NEST), base2(order,RING);
779  double maxang = 2.01*base.max_pixrad();
780  for (int m=0; m<nsamples/10; ++m)
781  {
782  I pix = I(rng.rand_uni()*base.Npix());
783  fix_arr<I,8> nb,nb2;
784  vec3 pixpt = base.pix2vec(pix);
785  base.neighbors(pix,nb);
786  base2.neighbors(base.nest2ring(pix),nb2);
787  for (int n=0; n<8; ++n)
788  if (nb[n]<0)
789  planck_assert(nb2[n]<0,"neighbor inconsistency");
790  else
791  planck_assert(base.nest2ring(nb[n])==nb2[n],"neighbor inconsistency");
792  sort(&nb[0],&nb[0]+8);
793  int check=0;
794  for (int n=0; n<8; ++n)
795  {
796  if (nb[n]<0)
797  ++check;
798  else
799  {
800  if (v_angle(base.pix2vec(nb[n]),pixpt)>maxang)
801  FAIL(cout<<" PROBLEM: order = "<<order<<", pix = "<<pix<<endl)
802  if ((n>0) && (nb[n]==nb[n-1]))
803  FAIL(cout<<" PROBLEM: order = "<<order<<", pix = "<<pix<<endl)
804  }
805  }
806  planck_assert((check<=1)||((order==0)&&(check<=2)),"too few neighbors");
807  }
808  }
809  for (I nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
810  {
811  T_Healpix_Base<I> base (nside,RING,SET_NSIDE);
812  double maxang = 2.01*base.max_pixrad();
813  for (int m=0; m<nsamples/10; ++m)
814  {
815  I pix = I(rng.rand_uni()*base.Npix());
816  fix_arr<I,8> nb;
817  vec3 pixpt = base.pix2vec(pix);
818  base.neighbors(pix,nb);
819  for (int n=0; n<8; ++n)
820  if ((nb[n]>=0) && (v_angle(base.pix2vec(nb[n]),pixpt)>maxang))
821  FAIL(cout<<" PROBLEM: nside = "<<nside<<", pix = "<<pix<<endl)
822  }
823  }
824  }
825 
826 void check_swap_scheme()
827  {
828  cout << "testing whether double swap_scheme() returns the original map"
829  << endl << "(for orders 0 to 10)." << endl;
830  for (int order=0; order<=10; ++order)
831  {
832  Healpix_Map<uint8> map(order,NEST);
833  for (int m=0; m<map.Npix(); ++m) map[m]=uint8(m&0xFF);
834  map.swap_scheme();
835  map.swap_scheme();
836  for (int m=0; m<map.Npix(); ++m)
837  if (map[m]!=(m&0xFF))
838  FAIL(cout<<" PROBLEM: order = "<<order<<", pix = "<<m<<endl)
839  }
840  }
841 
842 void check_issue_229 (Healpix_Ordering_Scheme scheme)
843  {
844  cout << "checking issue #229" << endl;
845  int order=8;
846  Healpix_Map<bool> map (order,scheme);
847  map.fill(false);
848  Healpix_Map<vec3> vmap(order,scheme);
849  for (int m=0; m<vmap.Npix(); ++m)
850  vmap[m]=vmap.pix2vec(m);
851  pointing ptg(halfpi-0.1,0);
852  double rad=0.1;
853  auto pixset = map.query_disc(ptg,rad);
854  vec3 vptg=ptg;
855  double cosrad=cos(rad);
856  for (tsize j=0; j<pixset.nranges(); ++j)
857  for (int i=pixset.ivbegin(j); i<pixset.ivend(j); ++i)
858  map[i] = true;
859  for (int i=0; i<map.Npix(); ++i)
860  {
861  bool inside = dotprod(vmap[i],vptg)>cosrad;
862  if (inside^map[i])
863  FAIL(cout << " PROBLEM: issue 229" << endl)
864  }
865  }
866 
867 void check_query_disc_strict (Healpix_Ordering_Scheme scheme)
868  {
869  cout << "testing whether all pixels found by query_disc() really" << endl
870  << "lie inside the disk (and vice versa)" << endl;
871  cout << "Ordering scheme: " << (scheme==RING ? "RING" : "NEST") << endl;
872  rng.seed(42);
873  for (int order=0; order<=5; ++order)
874  {
875  Healpix_Map<bool> map (order,scheme);
876  map.fill(false);
877  Healpix_Map<vec3> vmap(order,scheme);
878  for (int m=0; m<vmap.Npix(); ++m)
879  vmap[m]=vmap.pix2vec(m);
880  for (int m=0; m<100000; ++m)
881  {
882  pointing ptg;
883  random_dir (ptg);
884  double rad = pi/1 * rng.rand_uni();
885  auto pixset = map.query_disc(ptg,rad);
886  vec3 vptg=ptg;
887  double cosrad=cos(rad);
888  for (tsize j=0; j<pixset.nranges(); ++j)
889  for (int i=pixset.ivbegin(j); i<pixset.ivend(j); ++i)
890  map[i] = true;
891  for (int i=0; i<map.Npix(); ++i)
892  {
893  bool inside = dotprod(vmap[i],vptg)>cosrad;
894  if (inside^map[i])
895  FAIL(cout<<" PROBLEM: order = "<<order<<", ptg = "<<ptg<<endl)
896  }
897  for (tsize j=0; j<pixset.nranges(); ++j)
898  for (int i=pixset.ivbegin(j); i<pixset.ivend(j); ++i)
899  map[i] = false;
900  }
901  }
902  }
903 
904 template<typename I>void check_query_disc()
905  {
906  cout << "checking query_disc() " << bname<I>() << endl;
907  rng.seed(48);
908  int omax=min<int>(20,T_Healpix_Base<I>::order_max);
909  for (int order=0; order<=omax; ++order)
910  {
911  T_Healpix_Base<I> rbase (order,RING), nbase (order,NEST);
912  int niter=max(1,min(1000,100000>>order));
913  for (int m=0; m<niter; ++m)
914  {
915  pointing ptg;
916  random_dir (ptg);
917  double rad = pi/1 * rng.rand_uni();
918  auto pixset = rbase.query_disc(ptg,rad);
919  rangeset<I> pslast=pixset;
920  for (tsize fct=5; fct>0; --fct)
921  {
922  auto psi = rbase.query_disc_inclusive(ptg,rad,fct);
923  if (!psi.contains(pslast))
924  cout << " Potential problem: RING pixel sets inconsistent" << endl;
925  swap(pslast,psi);
926  }
927  I nval = pixset.nval();
928  pixset = nbase.query_disc(ptg,rad);
929  pslast=pixset;
930  for (tsize fct=8; fct>0; fct>>=1)
931  {
932  auto psi = nbase.query_disc_inclusive(ptg,rad,fct);
933  if (!psi.contains(pslast))
934  FAIL(cout << " PROBLEM: NEST pixel sets inconsistent" << endl)
935  swap(pslast,psi);
936  }
937  if (nval!=pixset.nval())
938  FAIL(cout << " PROBLEM: number of pixels different: "
939  << nval << " vs. " << pixset.nval() << endl)
940  }
941  }
942  }
943 template<typename I>void check_query_polygon()
944  {
945  cout << "checking query_polygon() " << bname<I>() << endl;
946  rng.seed(42);
947  int omax=min<int>(20,T_Healpix_Base<I>::order_max);
948  for (int order=0; order<=omax; ++order)
949  {
950  T_Healpix_Base<I> rbase (order,RING), nbase (order,NEST);
951  int niter=max(1,min(1000,100000>>order));
952  for (int m=0; m<niter; ++m)
953  {
954  vector<pointing> corner(3);
955  random_dir(corner[0]); random_dir(corner[1]); random_dir(corner[2]);
956  auto pixset = rbase.query_polygon(corner);
957  I nval = pixset.nval();
958  pixset = nbase.query_polygon(corner);
959  if (nval!=pixset.nval())
960  FAIL(cout << " PROBLEM: number of pixels different: "
961  << nval << " vs. " << pixset.nval() << endl)
962  pixset = rbase.query_polygon_inclusive(corner,4);
963  I nv1=pixset.nval();
964  pixset = nbase.query_polygon_inclusive(corner,4);
965  I nv2=pixset.nval();
966  if (nv1<nv2)
967  FAIL(cout << " PROBLEM: inclusive(RING)<inclusive(NEST): "
968  << nv1 << " vs. " << nv2 << endl)
969  if (nv2<nval)
970  FAIL(cout << " PROBLEM: inclusive(NEST)<non-inclusive: "
971  << nv2 << " vs. " << nval << endl)
972  }
973  }
974  }
975 
976 void helper_oop (int order)
977  {
978  Healpix_Map<double> map (order,RING), map2 (order,NEST), map3 (order,RING);
979  for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
980  map2.Import(map);
981  map3.Import(map2);
982  for (int m=0; m<map.Npix(); ++m)
983  if (!approx(map[m],map3[m],1e-12))
984  FAIL(cout << "PROBLEM: order = " << order << endl)
985  }
986 void helper_udgrade (int order, Healpix_Ordering_Scheme s1,
988  {
989  Healpix_Map<double> map (order,s1), map2 (order+2,s2), map3 (order, s1);
990  for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
991  map2.Import(map);
992  map3.Import(map2);
993  for (int m=0; m<map.Npix(); ++m)
994  if (!approx(map[m],map3[m],1e-12))
995  FAIL(cout << "PROBLEM: order = " << order << endl)
996  }
997 void helper_udgrade2 (int nside)
998  {
999  Healpix_Map<double> map (nside,RING,SET_NSIDE), map2 (nside*3,RING,SET_NSIDE),
1000  map3 (nside, RING,SET_NSIDE);
1001  for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
1002  map2.Import(map);
1003  map3.Import(map2);
1004  for (int m=0; m<map.Npix(); ++m)
1005  if (!approx(map[m],map3[m],1e-12))
1006  FAIL(cout << "PROBLEM: nside = " << nside << endl)
1007  }
1008 
1009 void check_import()
1010  {
1011  cout << "testing out-of-place swapping" << endl;
1012  for (int order=0; order<=7; ++order)
1013  helper_oop(order);
1014  cout << "testing downgrade(upgrade(map)) == map" << endl;
1015  for (int order=0; order<=7; ++order)
1016  {
1017  helper_udgrade(order,RING,RING);
1018  helper_udgrade(order,RING,NEST);
1019  helper_udgrade(order,NEST,NEST);
1020  helper_udgrade(order,NEST,RING);
1021  }
1022  for (int nside=3; nside<500; nside+=nside/2+1)
1023  helper_udgrade2(nside);
1024  }
1025 
1026 void check_average()
1027  {
1028  cout << "testing whether average(map) == average(downgraded map)" << endl;
1029  for (int order=1; order<=10; ++order)
1030  {
1031  Healpix_Map<double> map (order,RING), map2(1,RING);
1032  for (int m=0; m<map.Npix(); ++m)
1033  map[m] = rng.rand_uni()+0.01;
1034  map2.Import(map);
1035  double avg=map.average(), avg2=map2.average();
1036  if (!approx(avg,avg2,1e-13))
1037  FAIL(cout << "PROBLEM: order = " << order << " " << avg/avg2-1 << endl)
1038  }
1039  for (int nside=3; nside<1000; nside += nside/2+1)
1040  {
1041  Healpix_Map<double> map (nside,RING,SET_NSIDE), map2(1,RING,SET_NSIDE);
1042  for (int m=0; m<map.Npix(); ++m)
1043  map[m] = rng.rand_uni()+0.01;
1044  map2.Import(map);
1045  double avg=map.average(), avg2=map2.average();
1046  if (!approx(avg,avg2,1e-13))
1047  FAIL(cout << "PROBLEM: nside = " << nside << " " << avg/avg2-1 << endl)
1048  }
1049  }
1050 
1051 void random_alm (Alm<xcomplex<double> >&almT, Alm<xcomplex<double> >&almG,
1052  Alm<xcomplex<double> >&almC, int lmax, int mmax)
1053  {
1054  almT.Set(lmax,mmax); almG.Set(lmax,mmax); almC.Set(lmax,mmax);
1055 
1056  for (int l=0; l<=lmax; ++l)
1057  {
1058  almT(l,0)=dcomplex(rng.rand_gauss(),0.);
1059  almG(l,0)=dcomplex(rng.rand_gauss(),0.);
1060  almC(l,0)=dcomplex(rng.rand_gauss(),0.);
1061  }
1062 
1063  for (int m=1; m<=mmax; ++m)
1064  for (int l=m; l<=lmax; ++l)
1065  {
1066  almT(l,m).real(rng.rand_gauss()); almT(l,m).imag(rng.rand_gauss());
1067  almG(l,m).real(rng.rand_gauss()); almG(l,m).imag(rng.rand_gauss());
1068  almC(l,m).real(rng.rand_gauss()); almC(l,m).imag(rng.rand_gauss());
1069  }
1070  almG(0,0)=almC(0,0)=almG(1,0)=almC(1,0)=almG(1,1)=almC(1,1)=0;
1071  }
1072 
1073 void random_alm (Alm<xcomplex<double> >&alm, int lmax, int mmax)
1074  {
1075  alm.Set(lmax,mmax);
1076 
1077  for (int l=0; l<=lmax; ++l)
1078  { alm(l,0)=dcomplex(rng.rand_gauss(),0.); }
1079 
1080  for (int m=1; m<=mmax; ++m)
1081  for (int l=m; l<=lmax; ++l)
1082  { alm(l,m).real(rng.rand_gauss()); alm(l,m).imag(rng.rand_gauss()); }
1083  }
1084 
1085 void check_alm (const Alm<xcomplex<double> >&oalm,
1086  const Alm<xcomplex<double> >&alm, double epsilon)
1087  {
1088  for (int m=0; m<=alm.Mmax(); ++m)
1089  for (int l=m; l<=alm.Lmax(); ++l)
1090  {
1091  if (!abs_approx(oalm(l,m).real(),alm(l,m).real(),epsilon))
1092  FAIL(cout << "Problemr " << l << " " << m << endl)
1093  if (!abs_approx(oalm(l,m).imag(),alm(l,m).imag(),epsilon))
1094  FAIL(cout << "Problemi " << l << " " << m << endl)
1095  }
1096  }
1097 
1098 void check_alm2map2alm (int lmax, int mmax, int nside)
1099  {
1100  cout << "testing whether a_lm->map->a_lm returns original a_lm" << endl;
1101  cout << "lmax=" << lmax <<", mmax=" << mmax << ", nside=" << nside << endl;
1102  const double epsilon = 1e-8;
1103  Alm<xcomplex<double> > oalmT(lmax,mmax),almT(lmax,mmax),
1104  oalmG(lmax,mmax),almG(lmax,mmax),oalmC(lmax,mmax),almC(lmax,mmax);
1105  Healpix_Map<double> mapT(nside,RING,SET_NSIDE), mapQ(nside,RING,SET_NSIDE),
1106  mapU(nside,RING,SET_NSIDE);
1107 
1108  random_alm(oalmT,oalmG,oalmC,lmax,mmax);
1109  alm2map(oalmT,mapT);
1110  map2alm_iter2(mapT,almT,1e-12,1e-12);
1111  check_alm (oalmT, almT, epsilon);
1112 
1113  alm2map_spin(oalmG,oalmC,mapQ,mapU,1);
1114  map2alm_spin_iter2(mapQ,mapU,almG,almC,1,1e-12,1e-12);
1115  check_alm (oalmG, almG, epsilon);
1116  check_alm (oalmC, almC, epsilon);
1117 
1118  alm2map_pol(oalmT,oalmG,oalmC,mapT,mapQ,mapU);
1119  map2alm_pol_iter2(mapT,mapQ,mapU,almT,almG,almC,1e-12,1e-12);
1120  check_alm (oalmT, almT, epsilon);
1121  check_alm (oalmG, almG, epsilon);
1122  check_alm (oalmC, almC, epsilon);
1123  }
1124 
1125 void check_smooth_alm ()
1126  {
1127  cout << "testing whether unsmooth(smooth(a_lm)) returns a_lm" << endl;
1128  const double epsilon = 1e-14;
1129  const double fwhm = 100.*arcmin2rad;
1130  const int lmax=300, mmax=300;
1131  Alm<xcomplex<double> > oalmT(lmax,mmax),almT(lmax,mmax),
1132  oalmG(lmax,mmax),almG(lmax,mmax),oalmC(lmax,mmax),almC(lmax,mmax);
1133 
1134  random_alm(oalmT,oalmG,oalmC,lmax,mmax);
1135 
1136  almT=oalmT; almG=oalmG; almC=oalmC;
1137  smoothWithGauss (almT, fwhm);
1138  smoothWithGauss (almT, -fwhm);
1139  check_alm (oalmT, almT, epsilon);
1140  almT=oalmT;
1141  smoothWithGauss (almT, almG, almC, fwhm);
1142  smoothWithGauss (almT, almG, almC, -fwhm);
1143  check_alm (oalmT, almT, epsilon);
1144  check_alm (oalmG, almG, epsilon);
1145  check_alm (oalmC, almC, epsilon);
1146  }
1147 
1148 void check_rot_alm ()
1149  {
1150  cout << "testing whether rot^-1(rot(a_lm)) returns a_lm" << endl;
1151  const double epsilon = 2e-13;
1152  const int lmax=300;
1153  Alm<xcomplex<double> > oalm(lmax,lmax),alm(lmax,lmax);
1154 
1155  random_alm(oalm,lmax,lmax);
1156 
1157  alm=oalm;
1158  rotate_alm (alm,3,4,5);
1159  rotate_alm (alm,-5,-4,-3);
1160  check_alm (oalm, alm, epsilon);
1161  }
1162 
1163 double map_residual (const Healpix_Map<double> &a, const Healpix_Map<double> &b)
1164  {
1165  planck_assert(a.conformable(b), "maps are not conformable");
1166  double res=0;
1167  for (int i=0; i<a.Npix(); ++i)
1168  {
1169  double diff=a[i]-b[i];
1170  res+=diff*diff;
1171  }
1172  return res;
1173  }
1174 
1175 void check_ringweights ()
1176  {
1177  cout << "testing the accuracy of ring weights" << endl;
1178  int nside=128;
1179  double dummy;
1180  auto rwgt0=get_ringweights(nside,3*nside,1e-6,3000,dummy);
1181  arr<double> rwgt(2*nside);
1182  for (tsize i=0; i<rwgt.size(); ++i) rwgt[i]=rwgt0[i]+1.;
1183  Alm<xcomplex<double> > alm; random_alm(alm,1.5*nside,0);
1184  Healpix_Map<double> omap(nside,RING,SET_NSIDE),map2(omap);
1185  alm2map(alm,omap);
1186  map2alm(omap,alm,arr<double>(2*nside,1.));
1187  alm2map(alm,map2);
1188  double ressimple = map_residual(omap,map2);
1189  map2alm(omap,alm,rwgt);
1190  alm2map(alm,map2);
1191  double reswgt = map_residual(omap,map2);
1192  double resquot=sqrt(reswgt/ressimple);
1193  if (resquot>1e-5)
1194  FAIL(cout << " PROBLEM with ringweights: resquot = " << resquot << endl)
1195  }
1196 
1197 void check_fullweights ()
1198  {
1199  cout << "testing the accuracy of pixel weights" << endl;
1200  int nside=128;
1201  double dummy;
1202  auto fwgt=get_fullweights(nside,3*nside,1e-6,3000,dummy);
1203  Alm<xcomplex<double> > alm; random_alm(alm,1.5*nside,1.5*nside);
1204  Healpix_Map<double> omap(nside,RING,SET_NSIDE),map2(omap);
1205  alm2map(alm,omap);
1206  map2alm(omap,alm,arr<double>(2*nside,1.));
1207  alm2map(alm,map2);
1208  double ressimple = map_residual(omap,map2);
1209  map2=omap;
1210  apply_fullweights(map2,fwgt);
1211  map2alm(map2,alm,arr<double>(2*nside,1.));
1212  alm2map(alm,map2);
1213  double reswgt = map_residual(omap,map2);
1214  double resquot=sqrt(reswgt/ressimple);
1215  if (resquot>1e-5)
1216  FAIL(cout << " PROBLEM with full weights: resquot = " << resquot << endl)
1217  }
1218 
1219 void check_isqrt()
1220  {
1221  cout << "testing whether isqrt() works reliably" << endl;
1222  uint64 val=uint64(0xF234)<<16, valsq=val*val;
1223  if (isqrt(valsq)!=val) FAIL(cout << "PROBLEM1" << endl)
1224  if (isqrt(valsq-1)!=val-1) FAIL(cout << "PROBLEM2" << endl)
1225  }
1226 
1227 void check_pix2ang_acc()
1228  {
1229  cout << "testing accuracy of pix2ang at the poles" << endl;
1230  for (int m=0; m<=29;++m)
1231  {
1232  Healpix_Base2 base(m,RING);
1233  if (base.pix2ang(1).theta==0.)
1234  FAIL(cout << "PROBLEM: order " << m << endl)
1235  }
1236  }
1237 
1238 const int nsteps=1000000;
1239 
1240 template<typename I>void perf_neighbors(const string &name,
1241  Healpix_Ordering_Scheme scheme)
1242  {
1243  tsize cnt=0;
1244  wallTimers.start(name);
1246  for (int order=0; order<=omax; ++order)
1247  {
1248  T_Healpix_Base<I> base (order,scheme);
1249  I dpix=max(base.Npix()/nsteps,I(1));
1250  fix_arr<I,8> nres;
1251  for (I pix=0; pix<base.Npix(); pix+=dpix)
1252  {
1253  base.neighbors(pix,nres);
1254  ++cnt;
1255  }
1256  }
1257  wallTimers.stop(name);
1258  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1259  }
1260 template<typename I>void perf_pix2ang(const string &name,
1261  Healpix_Ordering_Scheme scheme, double &dummy)
1262  {
1263  tsize cnt=0;
1264  wallTimers.start(name);
1266  for (int order=0; order<=omax; ++order)
1267  {
1268  T_Healpix_Base<I> base (order,scheme);
1269  I dpix=max(base.Npix()/nsteps,I(1));
1270  for (I pix=0; pix<base.Npix(); pix+=dpix)
1271  {
1272  pointing p(base.pix2ang(pix));
1273  dummy+=p.theta+p.phi;
1274  ++cnt;
1275  }
1276  }
1277  wallTimers.stop(name);
1278  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1279  }
1280 template<typename I>void perf_pix2vec(const string &name,
1281  Healpix_Ordering_Scheme scheme, double &dummy)
1282  {
1283  tsize cnt=0;
1284  wallTimers.start(name);
1286  for (int order=0; order<=omax; ++order)
1287  {
1288  T_Healpix_Base<I> base (order,scheme);
1289  I dpix=max(base.Npix()/nsteps,I(1));
1290  for (I pix=0; pix<base.Npix(); pix+=dpix)
1291  {
1292  vec3 v(base.pix2vec(pix));
1293  dummy+=v.x+v.y+v.z;
1294  ++cnt;
1295  }
1296  }
1297  wallTimers.stop(name);
1298  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1299  }
1300 template<typename I>void perf_pix2zphi(const string &name,
1301  Healpix_Ordering_Scheme scheme, double &dummy)
1302  {
1303  tsize cnt=0;
1304  wallTimers.start(name);
1306  for (int order=0; order<=omax; ++order)
1307  {
1308  T_Healpix_Base<I> base (order,scheme);
1309  I dpix=max(base.Npix()/nsteps,I(1));
1310  double z,phi;
1311  for (I pix=0; pix<base.Npix(); pix+=dpix)
1312  {
1313  base.pix2zphi(pix,z,phi);
1314  dummy+=z+phi;
1315  ++cnt;
1316  }
1317  }
1318  wallTimers.stop(name);
1319  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1320  }
1321 template<typename I>void perf_zphi2pix(const string &name,
1322  Healpix_Ordering_Scheme scheme, double &dummy)
1323  {
1324  tsize cnt=0;
1325  wallTimers.start(name);
1327  double dz=2./sqrt(nsteps);
1328  double dph=twopi/sqrt(nsteps);
1329  for (int order=0; order<=omax; ++order)
1330  {
1331  T_Healpix_Base<I> base (order,scheme);
1332  for (double z=-1; z<1; z+=dz)
1333  for (double phi=0; phi<twopi; phi+=dph)
1334  {
1335  dummy+=base.zphi2pix(z,phi);
1336  ++cnt;
1337  }
1338  }
1339  wallTimers.stop(name);
1340  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1341  }
1342 template<typename I>void perf_ang2pix(const string &name,
1343  Healpix_Ordering_Scheme scheme, double &dummy)
1344  {
1345  tsize cnt=0;
1346  wallTimers.start(name);
1348  double dth=pi/sqrt(nsteps);
1349  double dph=twopi/sqrt(nsteps);
1350  for (int order=0; order<=omax; ++order)
1351  {
1352  T_Healpix_Base<I> base (order,scheme);
1353  for (double theta=0; theta<pi; theta+=dth)
1354  for (double phi=0; phi<twopi; phi+=dph)
1355  {
1356  dummy+=base.ang2pix(pointing(theta+1e-15*phi,phi));
1357  ++cnt;
1358  }
1359  }
1360  wallTimers.stop(name);
1361  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1362  }
1363 template<typename I>void perf_ring2nest(const string &name,double &dummy)
1364  {
1365  tsize cnt=0;
1366  wallTimers.start(name);
1368  for (int order=0; order<=omax; ++order)
1369  {
1370  T_Healpix_Base<I> base (order,RING);
1371  I dpix=max(base.Npix()/nsteps,I(1));
1372  for (I pix=0; pix<base.Npix(); pix+=dpix)
1373  {
1374  dummy+=base.ring2nest(pix);
1375  ++cnt;
1376  }
1377  }
1378  wallTimers.stop(name);
1379  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1380  }
1381 template<typename I>void perf_nest2ring(const string &name,double &dummy)
1382  {
1383  tsize cnt=0;
1384  wallTimers.start(name);
1386  for (int order=0; order<=omax; ++order)
1387  {
1388  T_Healpix_Base<I> base (order,RING);
1389  I dpix=max(base.Npix()/nsteps,I(1));
1390  for (I pix=0; pix<base.Npix(); pix+=dpix)
1391  {
1392  dummy+=base.nest2ring(pix);
1393  ++cnt;
1394  }
1395  }
1396  wallTimers.stop(name);
1397  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1398  }
1399 template<typename I>void perf_peano2nest(const string &name,double &dummy)
1400  {
1401  tsize cnt=0;
1402  wallTimers.start(name);
1404  for (int order=0; order<=omax; ++order)
1405  {
1406  T_Healpix_Base<I> base (order,NEST);
1407  I dpix=max(base.Npix()/nsteps,I(1));
1408  for (I pix=0; pix<base.Npix(); pix+=dpix)
1409  {
1410  dummy+=base.peano2nest(pix);
1411  ++cnt;
1412  }
1413  }
1414  wallTimers.stop(name);
1415  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1416  }
1417 template<typename I>void perf_nest2peano(const string &name,double &dummy)
1418  {
1419  tsize cnt=0;
1420  wallTimers.start(name);
1422  for (int order=0; order<=omax; ++order)
1423  {
1424  T_Healpix_Base<I> base (order,NEST);
1425  I dpix=max(base.Npix()/nsteps,I(1));
1426  for (I pix=0; pix<base.Npix(); pix+=dpix)
1427  {
1428  dummy+=base.nest2peano(pix);
1429  ++cnt;
1430  }
1431  }
1432  wallTimers.stop(name);
1433  cout << name << ": " << cnt/wallTimers.acc(name)*1e-6 << "MOps/s" << endl;
1434  }
1435 template<typename I>void perf_query_disc(const string &name,
1436  Healpix_Ordering_Scheme scheme, double &dummy)
1437  {
1438  tsize cnt=0;
1439  T_Healpix_Base<I> base(1024,scheme,SET_NSIDE);
1440  wallTimers.start(name);
1441  for (int m=0; m<1000; ++m)
1442  {
1443  dummy += base.query_disc(vec3(1,0,0),halfpi/9).nranges();
1444  ++cnt;
1445  }
1446  wallTimers.stop(name);
1447  cout << name << ": " << cnt/wallTimers.acc(name)*1e-3 << "kOps/s" << endl;
1448  }
1449 template<typename I>void perf_query_triangle(const string &name,
1450  Healpix_Ordering_Scheme scheme, double &dummy)
1451  {
1452  tsize cnt=0;
1453  T_Healpix_Base<I> base(1024,scheme,SET_NSIDE);
1454  vector<pointing> corner;
1455  corner.push_back(vec3(1,0.01,0.01));
1456  corner.push_back(vec3(0.01,1,0.01));
1457  corner.push_back(vec3(0.01,0.01,1));
1458  wallTimers.start(name);
1459  for (int m=0; m<1000; ++m)
1460  {
1461  dummy += base.query_polygon(corner).nranges();
1462  ++cnt;
1463  }
1464  wallTimers.stop(name);
1465  cout << name << ": " << cnt/wallTimers.acc(name)*1e-3 << "kOps/s" << endl;
1466  }
1467 template<typename I>void perf_query_polygon(const string &name,
1468  Healpix_Ordering_Scheme scheme, double &dummy)
1469  {
1470  tsize cnt=0;
1471  T_Healpix_Base<I> base(1024,scheme,SET_NSIDE);
1472  vector<pointing> corner;
1473  corner.push_back(vec3(1,0.01,0.01));
1474  corner.push_back(vec3(1,1,-0.3));
1475  corner.push_back(vec3(0.01,1,0.01));
1476  corner.push_back(vec3(0.01,0.01,1));
1477  wallTimers.start(name);
1478  for (int m=0; m<1000; ++m)
1479  {
1480  dummy += base.query_polygon(corner).nranges();
1481  ++cnt;
1482  }
1483  wallTimers.stop(name);
1484  cout << name << ": " << cnt/wallTimers.acc(name)*1e-3 << "kOps/s" << endl;
1485  }
1486 
1487 void perftest()
1488  {
1489  double dummy=0;
1490  cout << "Measuring performance of Healpix_Base methods." << endl;
1491  perf_pix2zphi<int> ("pix2zphi (RING):int ",RING,dummy);
1492  perf_pix2zphi<int> ("pix2zphi (NEST):int ",NEST,dummy);
1493  perf_pix2zphi<int64> ("pix2zphi (RING):int64",RING,dummy);
1494  perf_pix2zphi<int64> ("pix2zphi (NEST):int64",NEST,dummy);
1495  perf_zphi2pix<int> ("zphi2pix (RING):int ",RING,dummy);
1496  perf_zphi2pix<int> ("zphi2pix (NEST):int ",NEST,dummy);
1497  perf_zphi2pix<int64> ("zphi2pix (RING):int64",RING,dummy);
1498  perf_zphi2pix<int64> ("zphi2pix (NEST):int64",NEST,dummy);
1499  perf_pix2ang<int> ("pix2ang (RING):int ",RING,dummy);
1500  perf_pix2ang<int> ("pix2ang (NEST):int ",NEST,dummy);
1501  perf_pix2ang<int64> ("pix2ang (RING):int64",RING,dummy);
1502  perf_pix2ang<int64> ("pix2ang (NEST):int64",NEST,dummy);
1503  perf_ang2pix<int> ("ang2pix (RING):int ",RING,dummy);
1504  perf_ang2pix<int> ("ang2pix (NEST):int ",NEST,dummy);
1505  perf_ang2pix<int64> ("ang2pix (RING):int64",RING,dummy);
1506  perf_ang2pix<int64> ("ang2pix (NEST):int64",NEST,dummy);
1507  perf_pix2vec<int> ("pix2vec (RING):int ",RING,dummy);
1508  perf_pix2vec<int> ("pix2vec (NEST):int ",NEST,dummy);
1509  perf_pix2vec<int64> ("pix2vec (RING):int64",RING,dummy);
1510  perf_pix2vec<int64> ("pix2vec (NEST):int64",NEST,dummy);
1511  perf_neighbors<int> ("neighbors(NEST):int ",NEST);
1512  perf_neighbors<int> ("neighbors(RING):int ",RING);
1513  perf_neighbors<int64>("neighbors(NEST):int64",NEST);
1514  perf_neighbors<int64>("neighbors(RING):int64",RING);
1515  perf_ring2nest<int> ("ring2nest :int ",dummy);
1516  perf_ring2nest<int64>("ring2nest :int64",dummy);
1517  perf_nest2ring<int> ("nest2ring :int ",dummy);
1518  perf_nest2ring<int64>("nest2ring :int64",dummy);
1519  perf_peano2nest<int> ("peano2nest :int ",dummy);
1520  perf_peano2nest<int64>("peano2nest :int64",dummy);
1521  perf_nest2peano<int> ("nest2peano :int ",dummy);
1522  perf_nest2peano<int64>("nest2peano :int64",dummy);
1523  perf_query_disc<int> ("query_disc (RING):int ",RING,dummy);
1524  perf_query_disc<int> ("query_disc (NEST):int ",NEST,dummy);
1525  perf_query_disc<int64> ("query_disc (RING):int64",RING,dummy);
1526  perf_query_disc<int64> ("query_disc (NEST):int64",NEST,dummy);
1527  perf_query_triangle<int> ("query_triangle(RING):int ",RING,dummy);
1528  perf_query_triangle<int> ("query_triangle(NEST):int ",NEST,dummy);
1529  perf_query_triangle<int64>("query_triangle(RING):int64",RING,dummy);
1530  perf_query_triangle<int64>("query_triangle(NEST):int64",NEST,dummy);
1531  perf_query_polygon<int> ("query_polygon (RING):int ",RING,dummy);
1532  perf_query_polygon<int> ("query_polygon (NEST):int ",NEST,dummy);
1533  perf_query_polygon<int64> ("query_polygon (RING):int64",RING,dummy);
1534  perf_query_polygon<int64> ("query_polygon (NEST):int64",NEST,dummy);
1535 
1536  if (dummy<0) cout << dummy << endl;
1537  }
1538 
1539 } // unnamed namespace
1540 
1541 int main(int argc, const char **argv)
1542  {
1543  module_startup ("hpxtest",argc,argv,1,"");
1544  perftest();
1545  check_compress<int>();
1546  check_compress<unsigned>();
1547  check_compress<int64>();
1548  check_compress<uint64>();
1549  check_Moc<int>();
1550  check_Moc<int64>();
1551  check_rangeset<int>();
1552  check_rangeset<int64>();
1553  check_crangeset<int>();
1554  check_crangeset<int64>();
1555  check_isqrt();
1556  check_pix2ang_acc();
1557  check_smooth_alm();
1558  check_rot_alm();
1559  check_alm2map2alm(620,620,256);
1560  check_alm2map2alm(620,2,256);
1561  check_average();
1562  check_import();
1563  check_ringnestring<int>();
1564  check_ringnestring<int64>();
1565  check_nestpeanonest<int>();
1566  check_nestpeanonest<int64>();
1567  check_pixzphipix<int>();
1568  check_pixzphipix<int64>();
1569  check_zphipixzphi<int>();
1570  check_zphipixzphi<int64>();
1571  check_pixangpix<int>();
1572  check_pixangpix<int64>();
1573  check_neighbors<int>();
1574  check_neighbors<int64>();
1575  check_swap_scheme();
1576  check_query_disc_strict(RING);
1577  check_query_disc_strict(NEST);
1578  check_issue_229(RING);
1579  check_issue_229(NEST);
1580  check_query_disc<int>();
1581  check_query_disc<int64>();
1582  check_query_polygon<int>();
1583  check_query_polygon<int64>();
1584  check_ringweights();
1585  check_fullweights();
1586 #ifdef UNITTESTS
1587  if (errcount>0) planck_fail("unit test errors");
1588 #endif
1589  return 0;
1590  }
T nval() const
void append(const T &v1, const T &v2)
double v_angle(const vec3 &v1, const vec3 &v2)
vector< double > get_fullweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
bool overlaps(T a, T b) const
double rand_gauss()
unsigned int int_rand_uni()
bool contains(T a, T b) const
T nval() const
Definition: alm.h:88
rangeset op_xor(const rangeset &other) const
rangeset op_and(const rangeset &other) const
Healpix_Ordering_Scheme
void map2alm(const Healpix_Map< T > &map, Alm< xcomplex< T > > &alm, const arr< double > &weight, bool add_alm)
void apply_fullweights(Healpix_Map< T > &map, const std::vector< double > &wgt)
double phi
tsize nranges() const
void alm2map_pol(const Alm< xcomplex< T > > &almT, const Alm< xcomplex< T > > &almG, const Alm< xcomplex< T > > &almC, Healpix_Map< T > &mapT, Healpix_Map< T > &mapQ, Healpix_Map< T > &mapU, bool add_map)
void intersect(const T &a, const T &b)
std::size_t tsize
void seed(unsigned int x1=123456789, unsigned int y1=362436069, unsigned int z1=521288629, unsigned int w1=88675123)
void module_startup(const std::string &name, bool argc_valid, const std::string &usage, bool verbose=true)
bool overlaps(T a, T b) const
vec3_t< float64 > vec3
T dotprod(const vec3_t< T > &v1, const vec3_t< T > &v2)
void append(const T &v1, const T &v2)
bool approx(F a, F b, F epsilon=1e-5)
void alm2map(const Alm< xcomplex< T > > &alm, Healpix_Map< T > &map, bool add_map)
double rand_uni()
void add(const T &v1, const T &v2)
void remove(const T &v1, const T &v2)
rangeset op_or(const rangeset &other) const
const T & ivbegin(tdiff i) const
bool contains(T a, T b) const
#define planck_assert(testval, msg)
#define planck_fail(msg)
bool conformable(const T_Healpix_Base &other) const
Definition: healpix_base.h:435
const T & ivend(tdiff i) const
double theta
vector< double > get_ringweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
uint32 isqrt(I arg)
bool abs_approx(F a, F b, F epsilon=1e-5)
double cosdist_zphi(double z1, double phi1, double z2, double phi2)
I Npix() const
Definition: healpix_base.h:429
rangeset op_andnot(const rangeset &other) const

Generated on Mon Dec 10 2018 10:24:22 for Healpix C++