Skip to content

GitLab

Commit 718634a3 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

more functionality

parent 714edde6
Hide whitespace changes
Inline Side-by-side
Showing with 156 additions and 124 deletions
Healpix_cxx/healpix_map.cc
View file @ 718634a3
......@@ -37,44 +37,44 @@ namespace mr {
namespace healpix {
// template<typename T> void Healpix_Map<T>::Import_degrade
// (const Healpix_Map<T> &orig, bool pessimistic)
// {
// MR_assert(nside_<orig.nside_,"Import_degrade: this is no degrade");
// int fact = orig.nside_/nside_;
// MR_assert (orig.nside_==nside_*fact,
// "the larger Nside must be a multiple of the smaller one");
//
// int minhits = pessimistic ? fact*fact : 1;
// #pragma omp parallel
// {
// int m;
// #pragma omp for schedule (static)
// for (m=0; m<npix_; ++m)
// {
// int x,y,f;
// pix2xyf(m,x,y,f);
// int hits = 0;
// kahan_adder<double> adder;
// for (int j=fact*y; j<fact*(y+1); ++j)
// for (int i=fact*x; i<fact*(x+1); ++i)
// {
// int opix = orig.xyf2pix(i,j,f);
// if (!approx<double>(orig.map[opix],Healpix_undef))
// {
// ++hits;
// adder.add(orig.map[opix]);
// }
// }
// map[m] = T((hits<minhits) ? Healpix_undef : adder.result()/hits);
// }
// }
// }
//
// template void Healpix_Map<float>::Import_degrade
// (const Healpix_Map<float> &orig, bool pessimistic);
// template void Healpix_Map<double>::Import_degrade
// (const Healpix_Map<double> &orig, bool pessimistic);
template<typename T> void Healpix_Map<T>::Import_degrade
(const Healpix_Map<T> &orig, bool pessimistic)
{
MR_assert(nside_<orig.nside_,"Import_degrade: this is no degrade");
int fact = orig.nside_/nside_;
MR_assert (orig.nside_==nside_*fact,
"the larger Nside must be a multiple of the smaller one");
int minhits = pessimistic ? fact*fact : 1;
#pragma omp parallel
{
int m;
#pragma omp for schedule (static)
for (m=0; m<npix_; ++m)
{
int x,y,f;
pix2xyf(m,x,y,f);
int hits = 0;
tree_adder<double> adder;
for (int j=fact*y; j<fact*(y+1); ++j)
for (int i=fact*x; i<fact*(x+1); ++i)
{
int opix = orig.xyf2pix(i,j,f);
if (!approx<double>(orig.map[opix],Healpix_undef))
{
++hits;
adder.add(orig.map[opix]);
}
}
map[m] = T((hits<minhits) ? Healpix_undef : adder.result()/hits);
}
}
}
template void Healpix_Map<float>::Import_degrade
(const Healpix_Map<float> &orig, bool pessimistic);
template void Healpix_Map<double>::Import_degrade
(const Healpix_Map<double> &orig, bool pessimistic);
template<typename T> void Healpix_Map<T>::minmax (T &Min, T &Max) const
{
......
Healpix_cxx/healpix_map.h
View file @ 718634a3
......@@ -240,15 +240,15 @@ template<typename T> class Healpix_Map: public Healpix_Base
}
/*! Returns the average of all defined map pixels. */
// double average() const
// {
// kahan_adder<double> adder;
// int pix=0;
// for (int m=0; m<npix_; ++m)
// if (!approx<double>(map[m],Healpix_undef))
// { ++pix; adder.add(map[m]); }
// return (pix>0) ? adder.result()/pix : Healpix_undef;
// }
double average() const
{
tree_adder<double> adder;
int pix=0;
for (int m=0; m<npix_; ++m)
if (!approx<double>(map[m],Healpix_undef))
{ ++pix; adder.add(map[m]); }
return (pix>0) ? adder.result()/pix : Healpix_undef;
}
/*! Adds \a val to all defined map pixels. */
void Add (T val)
......
mr_util/math_utils.h
View file @ 718634a3
......@@ -34,6 +34,7 @@
#include <cmath>
#include <cstdint>
#include <vector>
namespace mr {
......@@ -175,6 +176,37 @@ template<typename T1, typename T2, typename... Args>
inline bool multiequal (const T1 &a, const T2 &b, Args... args)
{ return (a==b) && multiequal (a, args...); }
template<typename T> class tree_adder
{
private:
std::vector<T> state;
size_t n;
public:
tree_adder(): state(1,T(0)), n(0) {}
void add (const T &val)
{
state[0]+=val; ++n;
if (n>(size_t(1)<<(state.size()-1)))
state.push_back(T(0));
int shift=0;
while (((n>>shift)&1)==0)
{
state[shift+1]+=state[shift];
state[shift]=T(0);
++shift;
}
}
T result() const
{
T sum(0);
for (size_t i=0; i<state.size(); ++i)
sum+=state[i];
return sum;
}
};
/*! Returns \a atan2(y,x) if \a x!=0 or \a y!=0; else returns 0. */
inline double safe_atan2 (double y, double x)
{ return ((x==0.) && (y==0.)) ? 0.0 : std::atan2(y,x); }
......
test/hpxtest.cc
View file @ 718634a3
......@@ -1020,81 +1020,81 @@ template<typename I>void check_query_polygon()
}
}
// void helper_oop (int order)
// {
// Healpix_Map<double> map (order,RING), map2 (order,NEST), map3 (order,RING);
// for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
// map2.Import(map);
// map3.Import(map2);
// for (int m=0; m<map.Npix(); ++m)
// if (!approx(map[m],map3[m],1e-12))
// FAIL(cout << "PROBLEM: order = " << order << endl)
// }
// void helper_udgrade (int order, Ordering_Scheme s1,
// Ordering_Scheme s2)
// {
// Healpix_Map<double> map (order,s1), map2 (order+2,s2), map3 (order, s1);
// for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
// map2.Import(map);
// map3.Import(map2);
// for (int m=0; m<map.Npix(); ++m)
// if (!approx(map[m],map3[m],1e-12))
// FAIL(cout << "PROBLEM: order = " << order << endl)
// }
// void helper_udgrade2 (int nside)
// {
// Healpix_Map<double> map (nside,RING,SET_NSIDE), map2 (nside*3,RING,SET_NSIDE),
// map3 (nside, RING,SET_NSIDE);
// for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
// map2.Import(map);
// map3.Import(map2);
// for (int m=0; m<map.Npix(); ++m)
// if (!approx(map[m],map3[m],1e-12))
// FAIL(cout << "PROBLEM: nside = " << nside << endl)
// }
//
// void check_import()
// {
// cout << "testing out-of-place swapping" << endl;
// for (int order=0; order<=7; ++order)
// helper_oop(order);
// cout << "testing downgrade(upgrade(map)) == map" << endl;
// for (int order=0; order<=7; ++order)
// {
// helper_udgrade(order,RING,RING);
// helper_udgrade(order,RING,NEST);
// helper_udgrade(order,NEST,NEST);
// helper_udgrade(order,NEST,RING);
// }
// for (int nside=3; nside<500; nside+=nside/2+1)
// helper_udgrade2(nside);
// }
//
// void check_average()
// {
// cout << "testing whether average(map) == average(downgraded map)" << endl;
// for (int order=1; order<=10; ++order)
// {
// Healpix_Map<double> map (order,RING), map2(1,RING);
// for (int m=0; m<map.Npix(); ++m)
// map[m] = rng.rand_uni()+0.01;
// map2.Import(map);
// double avg=map.average(), avg2=map2.average();
// if (!approx(avg,avg2,1e-13))
// FAIL(cout << "PROBLEM: order = " << order << " " << avg/avg2-1 << endl)
// }
// for (int nside=3; nside<1000; nside += nside/2+1)
// {
// Healpix_Map<double> map (nside,RING,SET_NSIDE), map2(1,RING,SET_NSIDE);
// for (int m=0; m<map.Npix(); ++m)
// map[m] = rng.rand_uni()+0.01;
// map2.Import(map);
// double avg=map.average(), avg2=map2.average();
// if (!approx(avg,avg2,1e-13))
// FAIL(cout << "PROBLEM: nside = " << nside << " " << avg/avg2-1 << endl)
// }
// }
//
void helper_oop (int order)
{
Healpix_Map<double> map (order,RING), map2 (order,NEST), map3 (order,RING);
for (int m=0; m<map.Npix(); ++m) map[m] = frand()+0.01;
map2.Import(map);
map3.Import(map2);
for (int m=0; m<map.Npix(); ++m)
if (!approx(map[m],map3[m],1e-12))
FAIL(cout << "PROBLEM: order = " << order << endl)
}
void helper_udgrade (int order, Ordering_Scheme s1,
Ordering_Scheme s2)
{
Healpix_Map<double> map (order,s1), map2 (order+2,s2), map3 (order, s1);
for (int m=0; m<map.Npix(); ++m) map[m] = frand()+0.01;
map2.Import(map);
map3.Import(map2);
for (int m=0; m<map.Npix(); ++m)
if (!approx(map[m],map3[m],1e-12))
FAIL(cout << "PROBLEM: order = " << order << endl)
}
void helper_udgrade2 (int nside)
{
Healpix_Map<double> map (nside,RING,SET_NSIDE), map2 (nside*3,RING,SET_NSIDE),
map3 (nside, RING,SET_NSIDE);
for (int m=0; m<map.Npix(); ++m) map[m] = frand()+0.01;
map2.Import(map);
map3.Import(map2);
for (int m=0; m<map.Npix(); ++m)
if (!approx(map[m],map3[m],1e-12))
FAIL(cout << "PROBLEM: nside = " << nside << endl)
}
void check_import()
{
cout << "testing out-of-place swapping" << endl;
for (int order=0; order<=7; ++order)
helper_oop(order);
cout << "testing downgrade(upgrade(map)) == map" << endl;
for (int order=0; order<=7; ++order)
{
helper_udgrade(order,RING,RING);
helper_udgrade(order,RING,NEST);
helper_udgrade(order,NEST,NEST);
helper_udgrade(order,NEST,RING);
}
for (int nside=3; nside<500; nside+=nside/2+1)
helper_udgrade2(nside);
}
void check_average()
{
cout << "testing whether average(map) == average(downgraded map)" << endl;
for (int order=1; order<=10; ++order)
{
Healpix_Map<double> map (order,RING), map2(1,RING);
for (int m=0; m<map.Npix(); ++m)
map[m] = frand()+0.01;
map2.Import(map);
double avg=map.average(), avg2=map2.average();
if (!approx(avg,avg2,1e-13))
FAIL(cout << "PROBLEM: order = " << order << " " << avg/avg2-1 << endl)
}
for (int nside=3; nside<1000; nside += nside/2+1)
{
Healpix_Map<double> map (nside,RING,SET_NSIDE), map2(1,RING,SET_NSIDE);
for (int m=0; m<map.Npix(); ++m)
map[m] = frand()+0.01;
map2.Import(map);
double avg=map.average(), avg2=map2.average();
if (!approx(avg,avg2,1e-13))
FAIL(cout << "PROBLEM: nside = " << nside << " " << avg/avg2-1 << endl)
}
}
// void random_alm (Alm<xcomplex<double> >&almT, Alm<xcomplex<double> >&almG,
// Alm<xcomplex<double> >&almC, int lmax, int mmax)
// {
......@@ -1608,8 +1608,8 @@ int main(int argc, const char **argv)
// check_rot_alm();
// check_alm2map2alm(620,620,256);
// check_alm2map2alm(620,2,256);
// check_average();
// check_import();
check_average();
check_import();
check_ringnestring<int>();
check_ringnestring<int64_t>();
check_nestpeanonest<int>();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment