Healpix C++  3.50
needlet_tool_module.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) 2017 Max-Planck-Society
29  * Author: Martin Reinecke
30  */
31 
32 #include <memory>
33 #include <fstream>
34 #include "xcomplex.h"
35 #include "paramfile.h"
36 #include "alm.h"
37 #include "alm_fitsio.h"
38 #include "fitshandle.h"
39 #include "lsconstants.h"
40 #include "announce.h"
41 
42 using namespace std;
43 
44 namespace {
45 
46 class needlet_base
47  {
48  public:
49  virtual int lmin (int band) const = 0;
50  virtual int lmax (int band) const = 0;
51  virtual std::vector<double> getBand(int band) const = 0;
52  };
53 
54 
55 class needlet: public needlet_base
56  {
57  private:
58  enum { npoints=1000 };
59  std::vector<double> val;
60  double B;
61 
62  double phi(double t) const
63  {
64  using namespace std;
65  double pos=1.-2*B/(B-1.)*(t-1./B);
66  double p2=(val.size()-1)*abs(pos);
67  size_t ip2=size_t(p2);
68  if (ip2>=val.size()-1) return (pos>0.) ? 1. : 0.;
69  double delta=p2-ip2;
70  double v=val[ip2] + delta*(val[ip2+1]-val[ip2]);
71  return (pos>0) ? 0.5+v : 0.5-v;
72  }
73  double b2 (double xi) const
74  { return phi(xi/B) - phi(xi); }
75 
76  public:
77  needlet (double B_) : val(npoints), B(B_)
78  {
79  using namespace std;
80  val[0]=0;
81  double vold=exp(-1.);
82  for (int i=1; i<npoints-1; ++i)
83  {
84  double p0=double(i)/(npoints-1);
85  double v=exp(-1./(1-p0*p0));
86  val[i]=val[i-1]+0.5*(vold+v);
87  vold=v;
88  }
89  val[npoints-1]=val[npoints-2]+0.5*vold;
90  double xnorm=.5/val[npoints-1];
91  for (int i=0; i<npoints; ++i)
92  val[i]*=xnorm;
93  }
94 
95  virtual int lmin (int band) const
96  { return int(std::pow(B,band-1)+1.+1e-10); }
97  virtual int lmax (int band) const
98  { return int(std::pow(B,band+1)-1e-10); }
99  virtual std::vector<double> getBand(int band) const
100  {
101  using namespace std;
102  double fct=pow(B,-band);
103  vector<double> res(lmax(band)+1);
104  for (int i=0; i<lmin(band); ++i) res[i] = 0;
105  for (int i=lmin(band); i<int(res.size()); ++i)
106  res[i]=sqrt(b2(i*fct));
107  return res;
108  }
109  };
110 
111 class sdw: public needlet_base
112  {
113  private:
114  enum { npoints=1000 };
115  std::vector<double> val;
116  double B;
117 
118  double f_s2dw(double k)
119  {
120  if ((k<=1./B)||(k>=1.)) return 0;
121  double t = (k - (1. / B)) * (2.0 * B / (B-1.)) - 1.;
122  return exp(-2.0 / (1.0 - t*t)) / k;
123  }
124  double phi (double t) const
125  {
126  using namespace std;
127  double p2= (t-1./B)*(val.size()-1)/(1.-1./B);
128  int ip2=int(p2);
129  if (ip2>=int(val.size()-1)) return 0;
130  if (p2<=0.) return 1;
131  double delta=p2-ip2;
132  return 1.-(val[ip2] + delta*(val[ip2+1]-val[ip2]));
133  }
134  double b2 (double xi) const
135  { return phi(xi/B) - phi(xi); }
136 
137  public:
138  sdw (double B_) : val(npoints), B(B_)
139  {
140  using namespace std;
141  val[0]=0;
142  double vold=0.;
143  for (int i=1; i<npoints-1; ++i)
144  {
145  double p0=1./B+double(i)/(npoints-1)*(1-1./B);
146  double v=f_s2dw(p0);
147  val[i]=val[i-1]+0.5*(v+vold);
148  vold=v;
149  }
150  val[npoints-1]=val[npoints-2]+0.5*vold;
151  double xnorm=1./val[npoints-1];
152  for (int i=0; i<npoints; ++i)
153  val[i]*=xnorm;
154  }
155 
156  virtual int lmin (int band) const
157  { return int(std::pow(B,band-1)+1.+1e-10); }
158  virtual int lmax (int band) const
159  { return int(std::pow(B,band+1)-1e-10); }
160  virtual std::vector<double> getBand(int band) const
161  {
162  using namespace std;
163  double fct=pow(B,-band);
164  vector<double> res(lmax(band)+1);
165  for (int i=0; i<lmin(band); ++i) res[i] = 0;
166  for (int i=lmin(band); i<int(res.size()); ++i)
167  res[i]=sqrt(b2(i*fct));
168  return res;
169  }
170  };
171 
172 class spline: public needlet_base
173  {
174  private:
175  double B;
176 
177  static double psi (double x)
178  {
179  if (x<=0.) return 1.;
180  if (x>=1.) return 0.;
181  return (x>=0.5) ? 2.*(1.-x)*(1.-x)*(1.-x)
182  : 12.*x*x*(.5*x-.5) + 1.;
183  }
184 
185  double h (double xi) const
186  { return psi(xi/B)-psi(xi); }
187 
188  public:
189  spline (double B_) : B(B_) {}
190 
191  virtual int lmin (int /*band*/) const
192  { return 1; }
193  virtual int lmax (int band) const
194  { return int(std::pow(B,band+1)-1e-10); }
195  virtual std::vector<double> getBand(int band) const
196  {
197  using namespace std;
198  double fct=pow(B,-band);
199  vector<double> res(lmax(band)+1);
200  for (int i=0; i<lmin(band); ++i) res[i] = 0;
201  for (int i=lmin(band); i<int(res.size()); ++i)
202  res[i]=sqrt(h(i*fct));
203  return res;
204  }
205  };
206 
207 class cosine: public needlet_base
208  {
209  private:
210  std::vector<int> llim;
211 
212  void checkBand(int band) const
213  {
214  planck_assert((band>=0) && (band<int(llim.size()-2)), "illegal band");
215  }
216 
217  public:
218  cosine (const std::vector<int> &llim_) : llim(llim_)
219  {
220  planck_assert(llim.size()>=3, "llim vector needs at least 3 entries");
221  if ((llim[0]==0) && (llim[1]==0))
222  llim[0]=-1;
223  for (tsize i=1; i<llim.size(); ++i)
224  planck_assert(llim[i]>llim[i-1], "llim not strictly increasing");
225  }
226 
227  virtual int lmin (int band) const
228  {
229  checkBand(band);
230  return llim[band]+1;
231  }
232  virtual int lmax (int band) const
233  {
234  checkBand(band);
235  return llim[band+2]-1;
236  }
237  virtual std::vector<double> getBand(int band) const
238  {
239  using namespace std;
240  vector<double> res(lmax(band)+1);
241  for (int i=0; i<lmin(band); ++i) res[i]=0.;
242  for (int i=lmin(band); i<=llim[band+1]; ++i)
243  res[i]=cos(halfpi*(llim[band+1]-i)/(llim[band+1]-llim[band]));
244  for (int i=llim[band+1]; i<=lmax(band); ++i)
245  res[i]=cos(halfpi*(i-llim[band+1])/(llim[band+2]-llim[band+1]));
246  return res;
247  }
248  };
249 
250 template<typename T> void needlet_tool (paramfile &params)
251  {
252  string infile = params.template find<string>("infile");
253  string outfile = params.template find<string>("outfile");
254  string mode = params.template find<string>("mode");
255  planck_assert ((mode=="split") || (mode=="assemble"),"invalid mode");
256  bool split = mode=="split";
257  unique_ptr<needlet_base> needgen;
258  string ntype = params.template find<string>("needlet_type");
259  if (ntype=="cosine")
260  {
261  vector<string> lps=tokenize(params.template find<string>("llim"),',');
262  vector<int> lpeak;
263  for (auto x:lps) lpeak.push_back(stringToData<int>(x));
264  needgen.reset(new cosine(lpeak));
265  }
266  else if (ntype=="needatool")
267  needgen.reset(new needlet(params.template find<double>("B")));
268  else if (ntype=="sdw")
269  needgen.reset(new sdw(params.template find<double>("B")));
270  else if (ntype=="spline")
271  needgen.reset(new spline(params.template find<double>("B")));
272  else
273  planck_fail("unknown needlet type");
274  int loband=params.template find<int>("minband");
275  int hiband=params.template find<int>("maxband")+1;
276 
277  bool polarisation = params.template find<bool>("polarisation");
278  if (!polarisation)
279  {
280  if (split)
281  {
282  string outfile_needlets = params.template find<string>
283  ("outfile_needlets","");
284  int nlmax, nmmax;
285  get_almsize(infile, nlmax, nmmax);
286  auto alm = read_Alm_from_fits<T>(infile,nlmax,nmmax);
287  for (int i=loband; i<hiband; ++i)
288  {
289  int lmax_t=min(nlmax,needgen->lmax(i)),
290  mmax_t=min(nmmax, lmax_t);
291  Alm<xcomplex<T>> atmp(lmax_t,mmax_t);
292  for (int m=0; m<=mmax_t; ++m)
293  for (int l=m; l<=lmax_t; ++l)
294  atmp(l,m) = alm(l,m);
295  atmp.ScaleL(needgen->getBand(i));
296  write_Alm_to_fits (outfile+intToString(i,3)+".fits",atmp,lmax_t,mmax_t,
297  planckType<T>());
298  if (outfile_needlets!="")
299  {
300  fitshandle out;
301  out.create (outfile_needlets+intToString(i,3)+".fits");
302  vector<fitscolumn> cols;
303  cols.push_back(fitscolumn("coeff","[none]",1,PLANCK_FLOAT64));
304  out.insert_bintab(cols);
305  out.write_column(1,needgen->getBand(i));
306  }
307  }
308  }
309  else
310  {
311  Alm<xcomplex<T> > alm;
312  for (int i=loband; i<hiband; ++i)
313  {
314  int nlmax, nmmax;
315  get_almsize(infile+intToString(i,3)+".fits", nlmax, nmmax);
316  auto atmp = read_Alm_from_fits<T>(infile+intToString(i,3)+".fits",
317  nlmax,nmmax);
318  atmp.ScaleL(needgen->getBand(i));
319  for (int m=0; m<=alm.Mmax(); ++m)
320  for (int l=m; l<=alm.Lmax(); ++l)
321  atmp(l,m) += alm(l,m);
322  atmp.swap(alm);
323  write_Alm_to_fits (outfile,alm,alm.Lmax(),alm.Mmax(), planckType<T>());
324  }
325  }
326  }
327  else
328  {
329  planck_fail("polarisation not yet supported");
330  }
331  }
332 
333 } // unnamed namespace
334 
335 int needlet_tool (int argc, const char **argv)
336  {
337  module_startup ("needlet_tool", argc, argv);
338  paramfile params (getParamsFromCmdline(argc,argv));
339 
340  bool dp = params.find<bool> ("double_precision",false);
341  dp ? needlet_tool<double>(params) : needlet_tool<float>(params);
342 
343  return 0;
344  }
void create(const std::string &fname)
void write_Alm_to_fits(fitshandle &out, const Alm< xcomplex< T > > &alms, int lmax, int mmax, PDT datatype)
Definition: alm_fitsio.cc:151
void tokenize(const std::string &inp, char delim, std::vector< std::string > &list)
std::string intToString(int64 x, tsize width)
Definition: alm.h:88
void insert_bintab(const std::vector< fitscolumn > &cols, const std::string &extname="xtension")
int Lmax() const
Definition: alm.h:65
int Mmax() const
Definition: alm.h:67
T find(const std::string &key) const
void write_column(int colnum, const arr< T > &data, int64 offset=0)
paramfile getParamsFromCmdline(int argc, const char **argv, bool verbose=true)
std::size_t tsize
void module_startup(const std::string &name, bool argc_valid, const std::string &usage, bool verbose=true)
void get_almsize(fitshandle &inp, int &lmax, int &mmax)
Definition: alm_fitsio.cc:42
#define planck_assert(testval, msg)
void split(const std::string &inp, std::vector< T > &list)
#define planck_fail(msg)

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