75 #define FAIL(a) {a; if (++errcount>100) planck_fail("unit test errors");} 80 const int nsamples = 1000000;
91 void random_zphi (
double &z,
double &phi)
97 template<
typename I>
string bname()
98 {
return string(
"(basetype: ")+type2typename<I>()+
")"; }
101 template<
typename I>
rangeset<I> randomRangeSet(
int num, I start,
int dist)
105 for (
int i=0; i<num; ++i)
114 template<
typename I> Moc<I> randomMoc(
int num, I start,
int dist)
117 I curval=start+(I(1)<<(2*Moc<I>::maxorder));
118 for (
int i=0; i<num; ++i)
122 moc.addPixelRange(Moc<I>::maxorder,v1,v2);
128 template<
typename I>
rangeset<I> makeRS(
const string &input)
130 istringstream is(input);
158 template<
typename I>
void check_rangeset()
160 cout <<
"testing rangeset " << bname<I>() << endl;
177 planck_fail(
"Should have raised an IllegalArgumentException");
179 catch (PlanckError E) {}
218 planck_assert(b==makeRS<I>(
"-2 0 1 11 15 21 30 41"),
"error");
255 b = makeRS<I>(
"0 11 20 31");
258 b = makeRS<I>(
"0 11 20 31");
261 b = makeRS<I>(
"0 11 20 31");
264 b = makeRS<I>(
"0 11 20 31");
267 b = makeRS<I>(
"0 11 20 31");
273 makeRS<I>(
"20 31 40 51").op_or
274 (makeRS<I>(
"1 11 45 56")),
"error");
277 (makeRS<I>(
"1 11 45 56")),
"error");
279 makeRS<I>(
"1 11 45 56").op_or
280 (makeRS<I>(
"")),
"error");
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");
290 makeRS<I>(
"20 31 40 51").op_and
291 (makeRS<I>(
"")),
"error");
294 (makeRS<I>(
"20 31 40 51")),
"error");
298 makeRS<I>(
"20 31 40 51").op_andnot
299 (makeRS<I>(
"1 11 45 56")),
"error");
301 makeRS<I>(
"").op_andnot
302 (makeRS<I>(
"1 11 45 56")),
"error");
304 makeRS<I>(
"1 11 45 56").op_andnot
305 (makeRS<I>(
"")),
"error");
368 const int niter = 1000;
369 for (
int iter=0; iter<niter; ++iter)
380 template<
typename I>
crangeset<I> randomCRangeSet(
int num, I start,
int dist)
384 for (
int i=0; i<num; ++i)
393 template<
typename I>
crangeset<I> makeCRS(
const string &input)
395 istringstream is(input);
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");
423 template<
typename I>
void check_crangeset()
425 cout <<
"testing crangeset " << bname<I>() << endl;
442 planck_fail(
"Should have raised an IllegalArgumentException");
444 catch (PlanckError E) {}
462 makeCRS<I>(
"20 31 40 51").op_or
463 (makeCRS<I>(
"1 11 45 56")),
"error");
466 (makeCRS<I>(
"1 11 45 56")),
"error");
468 makeCRS<I>(
"1 11 45 56").op_or
469 (makeCRS<I>(
"")),
"error");
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");
479 makeCRS<I>(
"20 31 40 51").op_and
480 (makeCRS<I>(
"")),
"error");
482 makeCRS<I>(
"").op_and
483 (makeCRS<I>(
"20 31 40 51")),
"error");
487 makeCRS<I>(
"20 31 40 51").op_andnot
488 (makeCRS<I>(
"1 11 45 56")),
"error");
490 makeCRS<I>(
"").op_andnot
491 (makeCRS<I>(
"1 11 45 56")),
"error");
493 makeCRS<I>(
"1 11 45 56").op_andnot
494 (makeCRS<I>(
"")),
"error");
557 const int niter = 1000;
558 for (
int iter=0; iter<niter; ++iter)
569 template<
typename I>
void check_Moc()
571 cout <<
"testing MOC " << bname<I>() << endl;
573 moc.addPixelRange(0,4,5);
574 moc.addPixelRange(0,6,7);
575 moc.addPixelRange(2,4,17);
576 moc.addPixelRange(10,3000000,3000001);
581 Moc<I> xtmp = moc.degradedToOrder(8,
false);
585 xtmp=moc.degradedToOrder(8,
true);
589 planck_assert(moc==Moc<I>::fromCompressed(moc.toCompressed()),
"error");
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));
604 Moc<I> full; full.addPixelRange(0,0,12);
606 for (
tsize iter=0; iter<niter; ++iter)
608 Moc<I> a = randomMoc<I>(1000, 0, 100);
614 write_Moc_to_fits(
"!healpixtestmoctmp",a);
615 planck_assert(a==read_Moc_from_fits<I>(
"healpixtestmoctmp"),
"FITS problem");
620 template<
typename I>
void check_compress()
622 cout <<
"testing interpolation coding " << bname<I>() << endl;
626 for (
tsize x=0; x<100; ++x)
632 const vector<I> &v=b.data();
634 interpol_encode(v.begin(),v.end(),obs);
635 vector<uint8> comp=obs.state();
637 ibitstream ibs(comp);
638 interpol_decode(v2,ibs);
643 template<
typename I>
void check_ringnestring()
645 cout <<
"testing ring2nest(nest2ring(m))==m " << bname<I>() << endl;
646 for (
int order=0; order<=T_Healpix_Base<I>::order_max; ++order)
649 for (
int m=0; m<nsamples; ++m)
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)
658 template<
typename I>
void check_nestpeanonest()
660 cout <<
"testing peano2nest(nest2peano(m))==m " << bname<I>() << endl;
661 for (
int order=0; order<=T_Healpix_Base<I>::order_max; ++order)
664 for (
int m=0; m<nsamples; ++m)
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)
673 template<
typename I>
void check_pixzphipix()
675 cout <<
"testing zphi2pix(pix2zphi(m))==m " << bname<I>() << endl;
677 for (
int order=0; order<=omax; ++order)
680 for (
int m=0; m<nsamples; ++m)
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)
692 for (I nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
695 for (
int m=0; m<nsamples; ++m)
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)
706 template<
typename I>
void check_zphipixzphi()
708 cout <<
"testing pix2zphi(zphi2pix(ptg)) approx zphi " << bname<I>() << endl;
710 for (
int order=0; order<=omax; ++order)
713 double mincos = min (cos(base1.max_pixrad()),0.999999999999999);
714 for (
int m=0; m<nsamples; ++m)
716 double z,phi,z2,phi2;
718 base1.pix2zphi(base1.zphi2pix(z,phi),z2,phi2);
720 FAIL(cout <<
" PROBLEM: order = " << order
721 <<
", zphi = " << z <<
", " << phi << endl)
722 base2.pix2zphi(base2.zphi2pix(z,phi),z2,phi2);
724 FAIL(cout <<
" PROBLEM: order = " << order
725 <<
", zphi = " << z <<
", " << phi << endl)
728 for (
int nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
731 double mincos = min (cos(base.max_pixrad()),0.999999999999999);
732 for (
int m=0; m<nsamples; ++m)
734 double z,phi,z2,phi2;
736 base.pix2zphi(base.zphi2pix(z,phi),z2,phi2);
738 FAIL(cout <<
" PROBLEM: nside = " << nside
739 <<
", zphi = " << z <<
", " << phi << endl)
744 template<
typename I>
void check_pixangpix()
746 cout <<
"testing ang2pix(pix2ang(m))==m " << bname<I>() << endl;
748 for (
int order=0; order<=omax; ++order)
751 for (
int m=0; m<nsamples; ++m)
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)
760 for (I nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
763 for (
int m=0; m<nsamples; ++m)
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)
772 template<
typename I>
void check_neighbors()
774 cout <<
"testing neighbor function " << bname<I>() << endl;
776 for (
int order=0; order<=omax; ++order)
779 double maxang = 2.01*base.max_pixrad();
780 for (
int m=0; m<nsamples/10; ++m)
782 I pix = I(rng.
rand_uni()*base.Npix());
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)
791 planck_assert(base.nest2ring(nb[n])==nb2[n],
"neighbor inconsistency");
792 sort(&nb[0],&nb[0]+8);
794 for (
int n=0; n<8; ++n)
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)
806 planck_assert((check<=1)||((order==0)&&(check<=2)),
"too few neighbors");
809 for (I nside=3; nside<(I(1)<<omax); nside+=nside/2+1)
812 double maxang = 2.01*base.max_pixrad();
813 for (
int m=0; m<nsamples/10; ++m)
815 I pix = I(rng.
rand_uni()*base.Npix());
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)
826 void check_swap_scheme()
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)
833 for (
int m=0; m<map.Npix(); ++m) map[m]=uint8(m&0xFF);
836 for (
int m=0; m<map.Npix(); ++m)
837 if (map[m]!=(m&0xFF))
838 FAIL(cout<<
" PROBLEM: order = "<<order<<
", pix = "<<m<<endl)
844 cout <<
"checking issue #229" << endl;
849 for (
int m=0; m<vmap.Npix(); ++m)
850 vmap[m]=vmap.pix2vec(m);
853 auto pixset = map.query_disc(ptg,rad);
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)
859 for (
int i=0; i<map.Npix(); ++i)
861 bool inside =
dotprod(vmap[i],vptg)>cosrad;
863 FAIL(cout <<
" PROBLEM: issue 229" << endl)
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;
873 for (
int order=0; order<=5; ++order)
878 for (
int m=0; m<vmap.Npix(); ++m)
879 vmap[m]=vmap.pix2vec(m);
880 for (
int m=0; m<100000; ++m)
885 auto pixset = map.query_disc(ptg,rad);
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)
891 for (
int i=0; i<map.Npix(); ++i)
893 bool inside =
dotprod(vmap[i],vptg)>cosrad;
895 FAIL(cout<<
" PROBLEM: order = "<<order<<
", ptg = "<<ptg<<endl)
897 for (
tsize j=0; j<pixset.nranges(); ++j)
898 for (
int i=pixset.ivbegin(j); i<pixset.ivend(j); ++i)
904 template<
typename I>
void check_query_disc()
906 cout <<
"checking query_disc() " << bname<I>() << endl;
909 for (
int order=0; order<=omax; ++order)
912 int niter=max(1,min(1000,100000>>order));
913 for (
int m=0; m<niter; ++m)
918 auto pixset = rbase.query_disc(ptg,rad);
920 for (
tsize fct=5; fct>0; --fct)
922 auto psi = rbase.query_disc_inclusive(ptg,rad,fct);
923 if (!psi.contains(pslast))
924 cout <<
" Potential problem: RING pixel sets inconsistent" << endl;
927 I nval = pixset.nval();
928 pixset = nbase.query_disc(ptg,rad);
930 for (
tsize fct=8; fct>0; fct>>=1)
932 auto psi = nbase.query_disc_inclusive(ptg,rad,fct);
933 if (!psi.contains(pslast))
934 FAIL(cout <<
" PROBLEM: NEST pixel sets inconsistent" << endl)
937 if (nval!=pixset.nval())
938 FAIL(cout <<
" PROBLEM: number of pixels different: " 939 << nval <<
" vs. " << pixset.nval() << endl)
943 template<
typename I>
void check_query_polygon()
945 cout <<
"checking query_polygon() " << bname<I>() << endl;
948 for (
int order=0; order<=omax; ++order)
951 int niter=max(1,min(1000,100000>>order));
952 for (
int m=0; m<niter; ++m)
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);
964 pixset = nbase.query_polygon_inclusive(corner,4);
967 FAIL(cout <<
" PROBLEM: inclusive(RING)<inclusive(NEST): " 968 << nv1 <<
" vs. " << nv2 << endl)
970 FAIL(cout <<
" PROBLEM: inclusive(NEST)<non-inclusive: " 971 << nv2 <<
" vs. " << nval << endl)
976 void helper_oop (
int order)
979 for (
int m=0; m<map.Npix(); ++m) map[m] = rng.
rand_uni()+0.01;
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)
990 for (
int m=0; m<map.Npix(); ++m) map[m] = rng.
rand_uni()+0.01;
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)
997 void helper_udgrade2 (
int nside)
1000 map3 (nside,
RING,SET_NSIDE);
1001 for (
int m=0; m<map.Npix(); ++m) map[m] = rng.
rand_uni()+0.01;
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)
1011 cout <<
"testing out-of-place swapping" << endl;
1012 for (
int order=0; order<=7; ++order)
1014 cout <<
"testing downgrade(upgrade(map)) == map" << endl;
1015 for (
int order=0; order<=7; ++order)
1022 for (
int nside=3; nside<500; nside+=nside/2+1)
1023 helper_udgrade2(nside);
1026 void check_average()
1028 cout <<
"testing whether average(map) == average(downgraded map)" << endl;
1029 for (
int order=1; order<=10; ++order)
1032 for (
int m=0; m<map.Npix(); ++m)
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)
1039 for (
int nside=3; nside<1000; nside += nside/2+1)
1042 for (
int m=0; m<map.Npix(); ++m)
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)
1051 void random_alm (
Alm<xcomplex<double> >&almT,
Alm<xcomplex<double> >&almG,
1052 Alm<xcomplex<double> >&almC,
int lmax,
int mmax)
1054 almT.Set(lmax,mmax); almG.Set(lmax,mmax); almC.Set(lmax,mmax);
1056 for (
int l=0; l<=lmax; ++l)
1063 for (
int m=1; m<=mmax; ++m)
1064 for (
int l=m; l<=lmax; ++l)
1070 almG(0,0)=almC(0,0)=almG(1,0)=almC(1,0)=almG(1,1)=almC(1,1)=0;
1073 void random_alm (
Alm<xcomplex<double> >&alm,
int lmax,
int mmax)
1077 for (
int l=0; l<=lmax; ++l)
1080 for (
int m=1; m<=mmax; ++m)
1081 for (
int l=m; l<=lmax; ++l)
1085 void check_alm (
const Alm<xcomplex<double> >&oalm,
1086 const Alm<xcomplex<double> >&alm,
double epsilon)
1088 for (
int m=0; m<=alm.Mmax(); ++m)
1089 for (
int l=m; l<=alm.Lmax(); ++l)
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)
1098 void check_alm2map2alm (
int lmax,
int mmax,
int nside)
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;
1104 oalmG(lmax,mmax),almG(lmax,mmax),oalmC(lmax,mmax),almC(lmax,mmax);
1106 mapU(nside,
RING,SET_NSIDE);
1108 random_alm(oalmT,oalmG,oalmC,lmax,mmax);
1110 map2alm_iter2(mapT,almT,1e-12,1e-12);
1111 check_alm (oalmT, almT, epsilon);
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);
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);
1125 void check_smooth_alm ()
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;
1132 oalmG(lmax,mmax),almG(lmax,mmax),oalmC(lmax,mmax),almC(lmax,mmax);
1134 random_alm(oalmT,oalmG,oalmC,lmax,mmax);
1136 almT=oalmT; almG=oalmG; almC=oalmC;
1137 smoothWithGauss (almT, fwhm);
1138 smoothWithGauss (almT, -fwhm);
1139 check_alm (oalmT, almT, epsilon);
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);
1148 void check_rot_alm ()
1150 cout <<
"testing whether rot^-1(rot(a_lm)) returns a_lm" << endl;
1151 const double epsilon = 2e-13;
1155 random_alm(oalm,lmax,lmax);
1158 rotate_alm (alm,3,4,5);
1159 rotate_alm (alm,-5,-4,-3);
1160 check_alm (oalm, alm, epsilon);
1167 for (
int i=0; i<a.
Npix(); ++i)
1169 double diff=a[i]-b[i];
1175 void check_ringweights ()
1177 cout <<
"testing the accuracy of ring weights" << endl;
1182 for (
tsize i=0; i<rwgt.size(); ++i) rwgt[i]=rwgt0[i]+1.;
1188 double ressimple = map_residual(omap,map2);
1191 double reswgt = map_residual(omap,map2);
1192 double resquot=sqrt(reswgt/ressimple);
1194 FAIL(cout <<
" PROBLEM with ringweights: resquot = " << resquot << endl)
1197 void check_fullweights ()
1199 cout <<
"testing the accuracy of pixel weights" << endl;
1208 double ressimple = map_residual(omap,map2);
1213 double reswgt = map_residual(omap,map2);
1214 double resquot=sqrt(reswgt/ressimple);
1216 FAIL(cout <<
" PROBLEM with full weights: resquot = " << resquot << endl)
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)
1227 void check_pix2ang_acc()
1229 cout <<
"testing accuracy of pix2ang at the poles" << endl;
1230 for (
int m=0; m<=29;++m)
1233 if (base.pix2ang(1).theta==0.)
1234 FAIL(cout <<
"PROBLEM: order " << m << endl)
1238 const int nsteps=1000000;
1240 template<
typename I>
void perf_neighbors(
const string &name,
1244 wallTimers.start(name);
1246 for (
int order=0; order<=omax; ++order)
1249 I dpix=max(base.Npix()/nsteps,I(1));
1251 for (I pix=0; pix<base.Npix(); pix+=dpix)
1253 base.neighbors(pix,nres);
1257 wallTimers.stop(name);
1258 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1260 template<
typename I>
void perf_pix2ang(
const string &name,
1264 wallTimers.start(name);
1266 for (
int order=0; order<=omax; ++order)
1269 I dpix=max(base.Npix()/nsteps,I(1));
1270 for (I pix=0; pix<base.Npix(); pix+=dpix)
1273 dummy+=p.
theta+p.phi;
1277 wallTimers.stop(name);
1278 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1280 template<
typename I>
void perf_pix2vec(
const string &name,
1284 wallTimers.start(name);
1286 for (
int order=0; order<=omax; ++order)
1289 I dpix=max(base.Npix()/nsteps,I(1));
1290 for (I pix=0; pix<base.Npix(); pix+=dpix)
1292 vec3 v(base.pix2vec(pix));
1297 wallTimers.stop(name);
1298 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1300 template<
typename I>
void perf_pix2zphi(
const string &name,
1304 wallTimers.start(name);
1306 for (
int order=0; order<=omax; ++order)
1309 I dpix=max(base.Npix()/nsteps,I(1));
1311 for (I pix=0; pix<base.Npix(); pix+=dpix)
1313 base.pix2zphi(pix,z,phi);
1318 wallTimers.stop(name);
1319 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1321 template<
typename I>
void perf_zphi2pix(
const string &name,
1325 wallTimers.start(name);
1327 double dz=2./sqrt(nsteps);
1328 double dph=twopi/sqrt(nsteps);
1329 for (
int order=0; order<=omax; ++order)
1332 for (
double z=-1; z<1; z+=dz)
1333 for (
double phi=0; phi<twopi; phi+=dph)
1335 dummy+=base.zphi2pix(z,phi);
1339 wallTimers.stop(name);
1340 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1342 template<
typename I>
void perf_ang2pix(
const string &name,
1346 wallTimers.start(name);
1348 double dth=pi/sqrt(nsteps);
1349 double dph=twopi/sqrt(nsteps);
1350 for (
int order=0; order<=omax; ++order)
1353 for (
double theta=0; theta<pi; theta+=dth)
1354 for (
double phi=0; phi<twopi; phi+=dph)
1356 dummy+=base.ang2pix(
pointing(theta+1e-15*phi,phi));
1360 wallTimers.stop(name);
1361 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1363 template<
typename I>
void perf_ring2nest(
const string &name,
double &dummy)
1366 wallTimers.start(name);
1368 for (
int order=0; order<=omax; ++order)
1371 I dpix=max(base.Npix()/nsteps,I(1));
1372 for (I pix=0; pix<base.Npix(); pix+=dpix)
1374 dummy+=base.ring2nest(pix);
1378 wallTimers.stop(name);
1379 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1381 template<
typename I>
void perf_nest2ring(
const string &name,
double &dummy)
1384 wallTimers.start(name);
1386 for (
int order=0; order<=omax; ++order)
1389 I dpix=max(base.Npix()/nsteps,I(1));
1390 for (I pix=0; pix<base.Npix(); pix+=dpix)
1392 dummy+=base.nest2ring(pix);
1396 wallTimers.stop(name);
1397 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1399 template<
typename I>
void perf_peano2nest(
const string &name,
double &dummy)
1402 wallTimers.start(name);
1404 for (
int order=0; order<=omax; ++order)
1407 I dpix=max(base.Npix()/nsteps,I(1));
1408 for (I pix=0; pix<base.Npix(); pix+=dpix)
1410 dummy+=base.peano2nest(pix);
1414 wallTimers.stop(name);
1415 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1417 template<
typename I>
void perf_nest2peano(
const string &name,
double &dummy)
1420 wallTimers.start(name);
1422 for (
int order=0; order<=omax; ++order)
1425 I dpix=max(base.Npix()/nsteps,I(1));
1426 for (I pix=0; pix<base.Npix(); pix+=dpix)
1428 dummy+=base.nest2peano(pix);
1432 wallTimers.stop(name);
1433 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-6 <<
"MOps/s" << endl;
1435 template<
typename I>
void perf_query_disc(
const string &name,
1440 wallTimers.start(name);
1441 for (
int m=0; m<1000; ++m)
1443 dummy += base.query_disc(
vec3(1,0,0),halfpi/9).nranges();
1446 wallTimers.stop(name);
1447 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-3 <<
"kOps/s" << endl;
1449 template<
typename I>
void perf_query_triangle(
const string &name,
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)
1461 dummy += base.query_polygon(corner).nranges();
1464 wallTimers.stop(name);
1465 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-3 <<
"kOps/s" << endl;
1467 template<
typename I>
void perf_query_polygon(
const string &name,
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)
1480 dummy += base.query_polygon(corner).nranges();
1483 wallTimers.stop(name);
1484 cout << name <<
": " << cnt/wallTimers.acc(name)*1e-3 <<
"kOps/s" << endl;
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);
1536 if (dummy<0) cout << dummy << endl;
1541 int main(
int argc,
const char **argv)
1545 check_compress<int>();
1546 check_compress<unsigned>();
1547 check_compress<int64>();
1548 check_compress<uint64>();
1551 check_rangeset<int>();
1552 check_rangeset<int64>();
1553 check_crangeset<int>();
1554 check_crangeset<int64>();
1556 check_pix2ang_acc();
1559 check_alm2map2alm(620,620,256);
1560 check_alm2map2alm(620,2,256);
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();
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
unsigned int int_rand_uni()
bool contains(T a, T b) const
rangeset op_xor(const rangeset &other) const
rangeset op_and(const rangeset &other) const
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)
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)
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
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)
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)
bool conformable(const T_Healpix_Base &other) const
const T & ivend(tdiff i) const
vector< double > get_ringweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
bool abs_approx(F a, F b, F epsilon=1e-5)
double cosdist_zphi(double z1, double phi1, double z2, double phi2)
rangeset op_andnot(const rangeset &other) const