Healpix C++  3.50
weight_utils.cc
Go to the documentation of this file.
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 /*! \file weight_utils.cc
28  *
29  * Functionality for computing ring weights and full map weights
30  *
31  * Helpful literature:
32  * Graef, Kunis, Potts: On the computation of nonnegative quadrature weights
33  * on the sphere
34  * (https://www-user.tu-chemnitz.de/~potts/paper/quadgewS2.pdf)
35  * Lambers: Minimum Norm Solutions of Underdetermined Systems
36  * (http://www.math.usm.edu/lambers/mat419/lecture15.pdf)
37  * Shewchuk: An Introduction to the Conjugate Gradient Method Without the
38  * Agonizing Pain
39  * (https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf)
40  *
41  * Copyright (C) 2016-2018 Max-Planck-Society
42  * \author Martin Reinecke
43  */
44 
45 #include <vector>
46 #include <numeric>
47 #include <algorithm>
48 #include "healpix_map.h"
49 #include "alm_healpix_tools.h"
50 #include "weight_utils.h"
51 #include "alm.h"
52 #include "sharp_cxx.h"
53 #include "lsconstants.h"
54 
55 using namespace std;
56 
57 /* General considerations concerning the full weights and their a_lm
58  coefficients:
59  The Healpix grid geometry has several symmetries:
60  - mirror symmetry with respect to phi=0
61  - mirror symmetry with respect to the equator
62  - invariance with respect to rotation around the z-axis by multiples of
63  pi/2
64  This reduces the number of pixels with distinct weights to roughly 1/16 of
65  the total number of map pixels.
66  The geometrical symmetries also reduce the number of a_lm coefficients
67  required to describe the weight map:
68  - the mirror symmetry in phi means that all a_lm are purely real-valued
69  - the mirror symmetry in theta means that all a_lm with odd l are zero
70  - the rotational symmetry means that all a_lm with m not a multiple of 4
71  vanish as well.
72 
73  These symmetries are used when storing a_lm coefficients and weight maps to
74  save memory. So far, the SHTs themselves are carried out with full maps and
75  a_lm sets, because there is no libsharp support for this kind of symmetries
76  yet. */
77 
78 namespace {
79 
80 // inner product of two vectors
81 double dprod(const vector<double> &a, const vector<double> &b)
82  { return inner_product(a.begin(),a.end(),b.begin(),0.); }
83 
84 tsize n_fullweights (int nside)
85  { return ((3*nside+1)*(nside+1))/4; }
86 
87 template <typename T> void apply_weight (T &pix, double w, bool setwgt)
88  {
89  if (setwgt)
90  pix=T(w);
91  else
92  if (!approx<double>(pix,Healpix_undef)) pix*=T(1.+w);
93  }
94 
95 template <typename T> void apply_fullweights (Healpix_Map<T> &map,
96  const vector<double> &wgt, bool setwgt)
97  {
98  planck_assert (map.Scheme()==RING, "bad map ordering scheme");
99  int nside=map.Nside();
100  planck_assert(wgt.size()==n_fullweights(nside),
101  "incorrect size of weight array");
102  int pix=0, vpix=0;
103  for (int i=0; i<2*nside; ++i)
104  {
105  bool shifted = (i<nside-1) || ((i+nside)&1);
106  int qpix=min(nside,i+1);
107  bool odd=qpix&1;
108  int wpix=((qpix+1)>>1) + ((odd||shifted) ? 0 : 1);
109  int psouth=map.Npix()-pix-(qpix<<2);
110  for (int j=0; j<(qpix<<2); ++j)
111  {
112  int j4=j%qpix;
113  int rpix=min(j4,qpix - (shifted ? 1:0) - j4);
114  apply_weight (map[pix+j],wgt[vpix+rpix],setwgt);
115  if (i!=2*nside-1) // everywhere except on equator
116  apply_weight(map[psouth+j],wgt[vpix+rpix],setwgt);
117  }
118  pix+=qpix<<2;
119  vpix+=wpix;
120  }
121  }
122 
123 vector<double> extract_fullweights (const Healpix_Map<double> &map)
124  {
125  planck_assert (map.Scheme()==RING, "bad map ordering scheme");
126  int nside=map.Nside();
127  vector<double> res; res.reserve(n_fullweights(nside));
128  int pix=0;
129  for (int i=0; i<2*nside; ++i)
130  {
131  bool shifted = (i<nside-1) || ((i+nside)&1);
132  int qpix=min(nside,i+1);
133  bool odd=qpix&1;
134  int wpix=((qpix+1)>>1) + ((odd||shifted) ? 0 : 1);
135  for (int j=0; j<wpix; ++j)
136  res.push_back(map[pix+j]);
137  pix+=qpix<<2;
138  }
139  return res;
140  }
141 
142 tsize n_weightalm (int lmax, int mmax)
143  {
144  int nmsteps=(mmax>>2)+1;
145  return nmsteps * (((lmax+2)>>1) - nmsteps+1);
146  }
147 vector<double> extract_weightalm (const Alm<xcomplex<double> > &alm)
148  {
149  vector<double> res; res.reserve(n_weightalm(alm.Lmax(),alm.Mmax()));
150  for (int m=0; m<=alm.Mmax(); m+=4)
151  for (int l=m; l<=alm.Lmax(); l+=2)
152  res.push_back(real(alm(l,m)) * ((m==0) ? 1 : sqrt(2.)));
153  return res;
154  }
155 void expand_weightalm (const vector<double> &calm, Alm<xcomplex<double> > &alm)
156  {
157  planck_assert(calm.size()==n_weightalm(alm.Lmax(), alm.Mmax()),
158  "incorrect size of weight array");
159  alm.SetToZero();
160  tsize idx=0;
161  for (int m=0; m<=alm.Mmax(); m+=4)
162  for (int l=m; l<=alm.Lmax(); l+=2)
163  alm(l,m) = calm[idx++] * ((m==0) ? 1 : sqrt(0.5));
164  }
165 
166 class STS_hpwgt
167  {
168  private:
169  int lmax, mmax, nside;
170 
171  public:
172  using vectype=vector<double>;
173  STS_hpwgt (int lmax_, int mmax_, int nside_)
174  : lmax(lmax_), mmax(mmax_), nside(nside_)
175  { planck_assert((lmax&1)==0,"lmax must be even"); }
176  vectype S (const vectype &x) const
177  {
178  Alm<xcomplex<double> > ta(lmax,mmax);
179  expand_weightalm(x,ta);
180  Healpix_Map<double> tm(nside,RING, SET_NSIDE);
181  alm2map(ta,tm);
182  return extract_fullweights(tm);
183  }
184  vectype ST (const vectype &x) const
185  {
186  Healpix_Map<double> tm(nside,RING, SET_NSIDE);
187  apply_fullweights(tm,x,true);
188  Alm<xcomplex<double> > ta(lmax,mmax);
189  alm2map_adjoint(tm,ta);
190  return extract_weightalm(ta);
191  }
192  vectype apply (const vectype &x) const
193  { return ST(S(x)); }
194  };
195 
196 class STS_hpring
197  {
198  private:
199  int lmax, nside;
200  sharp_cxxjob<double> job;
201 
202  public:
203  using vectype=vector<double>;
204  STS_hpring (int lmax_, int nside_)
205  : lmax(lmax_), nside(nside_)
206  {
207  planck_assert((lmax&1)==0,"lmax must be even");
208  int nring=2*nside;
209  vector<double> dbl0(nring,0.),theta(nring);
210  vector<int> int1(nring,1);
211  vector<ptrdiff_t> ofs(nring);
212  Healpix_Base base(nside, RING, SET_NSIDE);
213  for (int i=0; i<nring; ++i)
214  {
215  ofs[i]=i;
216  int idum1,idum2;
217  bool bdum;
218  base.get_ring_info2 (i+1, idum1, idum2, theta[i], bdum);
219  }
220  job.set_general_geometry (nring, int1.data(), ofs.data(),
221  int1.data(), dbl0.data(), theta.data(), dbl0.data());
222  job.set_triangular_alm_info (lmax, 0);
223  }
224  vectype S(const vectype &alm) const
225  {
226  planck_assert(int(alm.size())==lmax/2+1,"bad input size");
227  vectype res(2*nside);
228  vector<xcomplex<double> > alm2(2*alm.size()-1,0.);
229  for (tsize i=0; i<alm.size(); ++i)
230  alm2[2*i]=alm[i];
231  job.alm2map(&alm2[0],&res[0],false);
232  return res;
233  }
234  vectype ST(const vectype &map) const
235  {
236  planck_assert(int(map.size())==2*nside,"bad input size");
237  vector<xcomplex<double> > alm2(lmax+1,0.);
238  job.alm2map_adjoint(&map[0],&alm2[0], false);
239  vectype res(lmax/2+1);
240  for (tsize i=0; i<res.size(); ++i)
241  res[i]=real(alm2[2*i]);
242  return res;
243  }
244  vectype apply (const vectype &x) const
245  { return ST(S(x)); }
246  };
247 
248 vector<double> muladd (double fct, const vector<double> &a,
249  const vector<double> &b)
250  {
251  planck_assert(a.size()==b.size(),"types not conformable");
252  vector<double> res(b);
253  for (tsize i=0; i<a.size(); ++i)
254  res[i] += fct*a[i];
255  return res;
256  }
257 
258 // canned algorithm B2 from Shewchuk
259 template<typename M> double cg_solve (const M &A, typename M::vectype &x,
260  const typename M::vectype &b, double epsilon, int itmax)
261  {
262  typename M::vectype r=muladd(-1.,A.apply(x),b), d(r);
263  double delta0=dprod(r,r), deltanew=delta0;
264  cout << "res0: " << sqrt(delta0) << endl;
265  for (int iter=0; iter<itmax; ++iter)
266  {
267  auto q=A.apply(d);
268  double alpha = deltanew/dprod(d,q);
269  x=muladd(alpha,d,x);
270  if (iter%300==0) // get accurate residual
271  r=muladd(-1.,A.apply(x),b);
272  else
273  r=muladd(-alpha,q,r);
274  double deltaold=deltanew;
275  deltanew=dprod(r,r);
276  cout << "\rIteration " << iter
277  << ": residual=" << sqrt(deltanew/delta0)
278  << " " << flush;
279  if (deltanew<epsilon*epsilon*delta0) { cout << endl; break; } // convergence
280  double beta=deltanew/deltaold;
281  d=muladd(beta,d,r);
282  }
283  return sqrt(deltanew/delta0);
284  }
285 
286 } // unnamed namespace
287 
288 vector<double> get_fullweights(int nside, int lmax, double epsilon, int itmax,
289  double &epsilon_out)
290  {
291  planck_assert((lmax&1)==0,"lmax must be even");
292  STS_hpwgt mat(lmax, lmax, nside);
293  vector<double> x(n_weightalm(lmax,lmax),0.);
294  vector<double> rhs=mat.ST(vector<double> (n_fullweights(nside),-1.));
295  rhs[0]+=12*nside*nside/sqrt(4*pi);
296  epsilon_out=cg_solve(mat, x, rhs, epsilon, itmax);
297  return mat.S(x);
298  }
299 
300 template <typename T> void apply_fullweights (Healpix_Map<T> &map,
301  const vector<double> &wgt)
302  { apply_fullweights (map,wgt,false); }
303 
304 template void apply_fullweights (Healpix_Map<float> &map,
305  const vector<double> &wgt);
306 template void apply_fullweights (Healpix_Map<double> &map,
307  const vector<double> &wgt);
308 
309 vector<double> get_ringweights(int nside, int lmax, double epsilon, int itmax,
310  double &epsilon_out)
311  {
312  planck_assert((lmax&1)==0,"lmax must be even");
313  STS_hpring mat(lmax,nside);
314  vector<double> nir(2*nside), x(lmax/2+1,0.);
315  for (tsize ith=0; ith<nir.size(); ++ith)
316  nir[ith]=8*min<int>(ith+1,nside);
317  nir[2*nside-1]/=2;
318  auto b=mat.ST(nir);
319  for (tsize i=0; i<b.size(); ++i)
320  b[i]=-b[i];
321  b[0]+=12*nside*nside/sqrt(4*pi);
322  epsilon_out=cg_solve(mat, x, b, epsilon, itmax);
323  auto mtmp=mat.S(x);
324  for (tsize i=0; i<mtmp.size(); ++i)
325  mtmp[i]/=nir[i];
326  return mtmp;
327  }
I Nside() const
Definition: healpix_base.h:427
vector< double > get_fullweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
Definition: alm.h:88
void apply_fullweights(Healpix_Map< T > &map, const std::vector< double > &wgt)
std::size_t tsize
void alm2map_adjoint(const Healpix_Map< T > &map, Alm< xcomplex< T > > &alm, bool add_alm)
Healpix_Ordering_Scheme Scheme() const
Definition: healpix_base.h:431
void alm2map(const Alm< xcomplex< T > > &alm, Healpix_Map< T > &map, bool add_map)
#define planck_assert(testval, msg)
const double Healpix_undef
Healpix value representing "undefined".
Definition: healpix_map.h:39
vector< double > get_ringweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
I Npix() const
Definition: healpix_base.h:429

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