52 #include "sharp_cxx.h" 81 double dprod(
const vector<double> &a,
const vector<double> &b)
82 {
return inner_product(a.begin(),a.end(),b.begin(),0.); }
84 tsize n_fullweights (
int nside)
85 {
return ((3*nside+1)*(nside+1))/4; }
87 template <
typename T>
void apply_weight (T &pix,
double w,
bool setwgt)
96 const vector<double> &wgt,
bool setwgt)
99 int nside=map.
Nside();
101 "incorrect size of weight array");
103 for (
int i=0; i<2*nside; ++i)
105 bool shifted = (i<nside-1) || ((i+nside)&1);
106 int qpix=min(nside,i+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)
113 int rpix=min(j4,qpix - (shifted ? 1:0) - j4);
114 apply_weight (map[pix+j],wgt[vpix+rpix],setwgt);
116 apply_weight(map[psouth+j],wgt[vpix+rpix],setwgt);
126 int nside=map.
Nside();
127 vector<double> res; res.reserve(n_fullweights(nside));
129 for (
int i=0; i<2*nside; ++i)
131 bool shifted = (i<nside-1) || ((i+nside)&1);
132 int qpix=min(nside,i+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]);
142 tsize n_weightalm (
int lmax,
int mmax)
144 int nmsteps=(mmax>>2)+1;
145 return nmsteps * (((lmax+2)>>1) - nmsteps+1);
147 vector<double> extract_weightalm (
const Alm<xcomplex<double> > &alm)
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.)));
155 void expand_weightalm (
const vector<double> &calm,
Alm<xcomplex<double> > &alm)
157 planck_assert(calm.size()==n_weightalm(alm.Lmax(), alm.Mmax()),
158 "incorrect size of weight array");
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));
169 int lmax, mmax, nside;
172 using vectype=vector<double>;
173 STS_hpwgt (
int lmax_,
int mmax_,
int nside_)
174 : lmax(lmax_), mmax(mmax_), nside(nside_)
176 vectype S (
const vectype &x)
const 179 expand_weightalm(x,ta);
182 return extract_fullweights(tm);
184 vectype ST (
const vectype &x)
const 190 return extract_weightalm(ta);
192 vectype apply (
const vectype &x)
const 200 sharp_cxxjob<double> job;
203 using vectype=vector<double>;
204 STS_hpring (
int lmax_,
int nside_)
205 : lmax(lmax_), nside(nside_)
209 vector<double> dbl0(nring,0.),theta(nring);
210 vector<int> int1(nring,1);
211 vector<ptrdiff_t> ofs(nring);
213 for (
int i=0; i<nring; ++i)
218 base.get_ring_info2 (i+1, idum1, idum2, theta[i], bdum);
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);
224 vectype S(
const vectype &alm)
const 227 vectype res(2*nside);
228 vector<xcomplex<double> > alm2(2*alm.size()-1,0.);
229 for (
tsize i=0; i<alm.size(); ++i)
231 job.alm2map(&alm2[0],&res[0],
false);
234 vectype ST(
const vectype &map)
const 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]);
244 vectype apply (
const vectype &x)
const 248 vector<double> muladd (
double fct,
const vector<double> &a,
249 const vector<double> &b)
252 vector<double> res(b);
253 for (
tsize i=0; i<a.size(); ++i)
259 template<
typename M>
double cg_solve (
const M &A,
typename M::vectype &x,
260 const typename M::vectype &b,
double epsilon,
int itmax)
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)
268 double alpha = deltanew/dprod(d,q);
271 r=muladd(-1.,A.apply(x),b);
273 r=muladd(-alpha,q,r);
274 double deltaold=deltanew;
276 cout <<
"\rIteration " << iter
277 <<
": residual=" << sqrt(deltanew/delta0)
279 if (deltanew<epsilon*epsilon*delta0) { cout << endl;
break; }
280 double beta=deltanew/deltaold;
283 return sqrt(deltanew/delta0);
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);
301 const vector<double> &wgt)
305 const vector<double> &wgt);
307 const vector<double> &wgt);
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);
319 for (
tsize i=0; i<b.size(); ++i)
321 b[0]+=12*nside*nside/sqrt(4*pi);
322 epsilon_out=cg_solve(mat, x, b, epsilon, itmax);
324 for (
tsize i=0; i<mtmp.size(); ++i)
vector< double > get_fullweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
void apply_fullweights(Healpix_Map< T > &map, const std::vector< double > &wgt)
void alm2map_adjoint(const Healpix_Map< T > &map, Alm< xcomplex< T > > &alm, bool add_alm)
Healpix_Ordering_Scheme Scheme() const
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".
vector< double > get_ringweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)