Commit a47ade7a authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

nilc

parent c1fde18b
This diff is collapsed.
/*
* Project: SZ
* Date: 08/10
* Authors: M. Le Jeune
*
*/
/* Copyright (C) 2008 APC CNRS Université Paris Diderot <lejeune@apc.univ-paris7.fr>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see http://www.gnu.org/licenses/gpl.html
*/
#include <cmath>
#include <iostream>
#include <string>
//#include "qromb.h"
#include <sstream>
#include "arr.h"
#include "vec3.h"
#include <gsl/gsl_integration.h>
using namespace std;
/// \param theta_min : min bound for integration (max is pi/2)
//double theta_min;
/// \param beta : beta profile parameter
//double beta = (double)2/3;
/// \param rmax : viriel radius
//double rmax;
/// \param theta_c : core radius
//double theta_c;
/// beta function to be integrated between theta_min and pi/2
/// \param theta : angular distance from cluster center
double betat (double theta, void * params);
/// beta function at center to be integrated between 0 and rmax
/// \param theta : angular distance from cluster center
double beta0 (double theta, void * params);
/// give the beta profile value.
/// \param theta_i : angular distance from cluster center
/// \param rmax : viriel radius
/// \param r_core : core radius
double beta_profile (double rmax, double r_core, double theta_i,double beta);
/// nwf function to be integrated between theta_min and pi/2
/// \param theta : angular distance from cluster center
double nwft (double theta, void * params);
/// nwf function at center to be integrated between 0 and rmax
/// \param theta : angular distance from cluster center
double nwf0 (double theta, void * params);
double xmm_profile (double rmax, double scale_radius, double theta_i, double alpha, double beta, double gamma);
/// nwf function at center to be integrated between 0 and rmax
/// \param theta : angular distance from cluster center
double nwf5 (double theta, void * params);
double y5r500 (double scale_radius, double alpha, double beta, double gamma);
## ----------------------------------------------------------------------------
## Set your config here
LIBS = [ 'healpix_cxx' ,'cxxsupport' ,'cfitsio', 'spherelib','processoptions', 'gsl', 'gslcblas','libpsht' ,'c_utils' ,'fftpack']
SPHERELIB = '../lib/src'
SPHEREINC = '../include'
HEALPIXLIB = '../healpix/spherelib/lib'
HEALPIXINC = '../healpix/spherelib/include'
HEALPIXDATA = '\\\"../healpix/data\\\"'
CFITSIOLIB = '/usr/local/lib'
CFITSIOINC = '/usr/local/include'
CXX = 'g++'
CXXFLAGS = '-w -O3 -fPIC -ffast-math -fopenmp -g0 -s'
## ----------------------------------------------------------------------------
import glob
cfiles = glob.glob('./*.cc')
opt = Environment(CXXFLAGS = CXXFLAGS,LIBS=LIBS, CXX=CXX,CPPDEFINES={'HEALPIXDATA': HEALPIXDATA}, CPPPATH = [SPHERELIB, SPHEREINC, HEALPIXINC,CFITSIOINC],LINKFLAGS=CXXFLAGS, LIBPATH=[SPHERELIB,HEALPIXLIB,CFITSIOLIB,"."] )
## Build healpix
print('building healpix (this may take a while)')
import os
target = "spherelib"
cfitsio = "/usr/local"
compil_result = os.system('cd ../healpix && '
'HEALPIX_TARGET=%s CFITSIO_EXT_PREFIX=%s make '%(target,cfitsio))
## Build the library
SConscript(SPHERELIB+'/Sconscript', exports='opt')
opt.Library('spherelib', glob.glob('../lib/src/*.cpp'))
opt.Library('processoptions', ['processoptions.cpp', 'processoptions.h'])
## Compile binaries
for p in range(len(cfiles)):
opt.Install('/usr/local/bin',opt.Program(cfiles[p]) )
opt.Alias('install', '/usr/local/bin')
/*
* Project: nilc
* Date: 11/01
* Authors: M. Le Jeune
*
*/
/* Copyright (C) 2011 APC CNRS Université Paris Diderot <lejeune@apc.univ-paris7.fr>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see http://www.gnu.org/licenses/gpl.html
*/
#include "fits_tool.h"
#include <iostream>
#include <string>
#include "paramfile.h"
#include "processoptions.h"
#include <sys/time.h>
int main(int argc ,const char** argv)
{
timeval start, end;
char** ptr ;
ptr = const_cast <char **> (argv);
/// options
bool almopt=false; // input are alm files
bool pixel_based=false; // pixel domain localisation only
/// inputs
string outfile; // healpix map FITS file
string *infile; // healpix map FITS file
int nbin;
/// parameters
arr<double> mixcol; // size is nbin
arr<int> bands; // not used
string *zonefile; // healpix map FITS file
int nbands;
int lmax;
int nside = -1;
int n_iter = 0;
// command line parsing (1/2)
string longopt[] = {"almopt", "n_iter", "nside"};
ArgType type[] = {FLAG, INT, INT};
string usage[] = {"Inputs are alm FITS file.", "Number of iteration of anafast.", "Resolution of the output map."};
char shortopt[] = "ain";
Args args = Args(3, longopt, shortopt, type, usage,"Perform needlet ILC.", "outfile l0 J zonefile1 lmax1 ... zonefileJ lmaxJ N infile1 coeff1 ... infileN coeffN ");
// parameter file
if ((argc == 2) && (argv[1][0]!='-'))
{
module_startup ("nilc", argc, argv, 2, "<parameter file>");
paramfile params (argv[1]);
almopt = params.find<bool> ("almopt",false);
nside = params.find<int> ("nside" ,-1);
n_iter = params.find<int> ("n_iter",0);
outfile= params.find<string> ("outfile","");
string sinfile = params.find<string> ("infile","");
nbin = parse_list (sinfile, &infile);
string smixcol = params.find<string> ("mixcol","");
string*stmixcol;
int nbin2 = parse_list (smixcol, &stmixcol);
if (nbin!=nbin2){
fprintf(stderr,"\nERROR!!! Can't parse 'mixcol' parameter. Found %d elements vs %d inputs\n",nbin2,nbin);
return 1;
}
mixcol.alloc (nbin);
for (int i=0; i<nbin; i++)
mixcol[i] = atof(stmixcol[i].c_str());
delete []stmixcol;
string szonefile = params.find<string> ("zonefile","");
nbands = parse_list (szonefile, &zonefile);
string sbands = params.find<string> ("bands","");
string* stbands;
int nbands2 = parse_list (sbands, &stbands);
if (nbands2 == 0) {
pixel_based = true;
nbands = 1;
}
else if (nbands!=nbands2-2){
fprintf(stderr,"\nERROR!!! Can't parse 'bands' parameter. Found %d elements vs %d zones\n",nbands2,nbands);
return 1;
}
else {
bands.alloc (nbands2);
for (int i=0; i<nbands2; i++)
bands[i] = atoi(stbands[i].c_str());
delete []stbands;
}
}
// command line parsing (2/2)
else {
args.processargs(argc, ptr);
args.info(); //args.printusage();
almopt = args["almopt"];
n_iter = args["n_iter"];
nside = args["nside"];
int idx = 1;
outfile = argv[idx++];
int l0 = atoi(argv[idx++]);
nbands = atoi(argv[idx++]);
bands.alloc(nbands+2);
zonefile = new string[nbands];
bands[0] = 0;
bands[1] = l0;
for (int i=0; i<nbands; i++){
zonefile[i] = argv[idx++];
bands[i+2] = atoi(argv[idx++]);}
nbin = atoi(argv[idx++]);
infile = new string[nbin];
mixcol.alloc(nbin);
for (int i=0; i<nbin; i++){
infile[i] = argv[idx++];
mixcol[i] = atof(argv[idx++]);}
}
// preprocessing
fitshandle inp;
if (!almopt){
inp.open(infile[0]);
inp.goto_hdu(2);
nside = (int) sqrt(inp.nelems(1)/12);
}
if (nside<0){
fprintf(stderr,"\nERROR!!! Missing 'nside' parameter\n");
return 1;
}
if (nbands==0){
get_lst_tmp (1, "zone", &zonefile);
MapExt<int> mask;
mask.SetNside(nside, RING);
mask.fill(0);
fitshandle out;
out.create(zonefile[0]);
write_Healpix_map_to_fits (out, mask,planckType<int>());
}
arr<double> onearr(nbands);
onearr.fill(1);
arr<int> nsides(nbands);
for (int j=0; j<nbands; j++) {
inp.open(zonefile[j]);
inp.goto_hdu(2);
nsides[j] = (int) sqrt(inp.nelems(1)/12);
}
// pixel based case
if (pixel_based){
string statfile = get_tmp_file("stat");
string ilcfile = get_tmp_file("ilc");
fprintf(stdout, "Compute second order stats: ");
gettimeofday(&start, NULL);
if (almopt){
string* mapfile;
get_lst_tmp (nbin, "map", &mapfile);
for (int m=0; m<nbin; m++)
alm2map<double> (infile[m], mapfile[m], nside, false);
map2stat<double> (mapfile, nbin, statfile, zonefile[0]);}
else
map2stat<double> (infile, nbin, statfile, zonefile[0]);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
fprintf(stdout, "Compute ILC filter: ");
gettimeofday(&start, NULL);
stat2ILCfilter<double> (statfile, nbin, ilcfile, mixcol);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
fprintf(stdout, "Apply ILC filter: ");
gettimeofday(&start, NULL);
combinezone<double> (infile, nbin, outfile, ilcfile,zonefile[0]);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
}
else{
string *winfile;
get_lst_tmp (nbands, "window", &winfile);
string* scalefile;
get_lst_tmp (nbands, "cfalm", &scalefile);
makesplinefilters (bands, winfile, 0);
for (int j=0; j<nbands; j++) {
string* mapfile;
get_lst_tmp (nbin, "needlet", &mapfile);
for (int m=0; m<nbin; m++) {
string almfile1 = get_tmp_file("ialm");
string almfile2 = get_tmp_file("falm");
fprintf(stdout, "Perform analysis of map %d and scale %d. Lmax is %d, Nside is %d: ", m, j, bands[j+2], nsides[j]);
gettimeofday(&start, NULL);
if (almopt)
filteralm<double> (infile[m], almfile2, winfile[j], false);
else{
map2alm<double> (infile[m], almfile1, bands[j+2], n_iter, "", false);
filteralm<double> (almfile1, almfile2, winfile[j], false);
}
alm2map<double> (almfile2, mapfile[m], nsides[j], false);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
}
string statfile = get_tmp_file("stat");
string ilcfile = get_tmp_file("ilc");
string cmapfile = get_tmp_file("cmap");
string calmfile = get_tmp_file("calm");
fprintf(stdout, "Compute ILC filter for scale %d: ", j);
gettimeofday(&start, NULL);
map2stat<double> (mapfile, nbin, statfile, zonefile[j]);
stat2ILCfilter<double> (statfile, nbin, ilcfile, mixcol);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
fprintf(stdout, "Apply ILC filter for scale %d: ", j);
gettimeofday(&start, NULL);
combinezone<double> (mapfile, nbin, cmapfile, ilcfile,zonefile[j]);
map2alm<double> (cmapfile, calmfile, bands[j+2-1], n_iter, "", false);
filteralm<double> (calmfile, scalefile[j], winfile[j], false);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
}
string almfile = get_tmp_file("alm");
fprintf(stdout, "Combine scales: ");
gettimeofday(&start, NULL);
combinealm<double>(scalefile, nbands, almfile, onearr, false);
alm2map<double> (almfile, outfile, nside, 0);
gettimeofday(&end, NULL);
print_time (start.tv_sec, start.tv_usec, end.tv_sec, end.tv_usec);
}
return 0;
}
// Copyright (C) 2008 Marc Betoule
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
// report bugs at betoule@apc.univ-paris7.fr
#include "processoptions.h"
#include <iostream>
#include <sstream>
#include <cstdlib>
Args::Args(int noptions, string * longoptions, char * shortoptions, ArgType * type, string * usage, string _purpose, string _fullusage){
longOpts = new struct option [noptions+1];
shortOpts = "";
optionroles = "";
for (int i = 0; i < noptions; i++){
types[shortoptions[i]] = type[i];
longnames[shortoptions[i]] = longoptions[i];
shortnames[longoptions[i]] = shortoptions[i];
isset[longoptions[i]] = false;
switch (type[i]){
case STRING :
{
stringvalues[longoptions[i]] = "";
struct option o = {longoptions[i].c_str(), required_argument, NULL, shortoptions[i]};
longOpts[i] = o;
shortOpts += shortoptions[i];
shortOpts += ":";
optionroles += "\n--"+longoptions[i]+", -"+shortoptions[i]+ " <string>\t: " + usage[i];
}
break;
case FLOAT :
{
floatvalues[longoptions[i]] = 0.0;
struct option o = {longoptions[i].c_str(), required_argument, NULL, shortoptions[i]};
longOpts[i] = o;
shortOpts += shortoptions[i];
shortOpts += ":";
optionroles += "\n--"+longoptions[i]+", -"+shortoptions[i]+ " <float>\t: " + usage[i];
}
break;
case FLAG :
{
flagvalues[longoptions[i]] = 0;
struct option o = {longoptions[i].c_str(), no_argument, NULL, shortoptions[i]};
longOpts[i] = o;
shortOpts += shortoptions[i];
optionroles += "\n--"+longoptions[i]+", -"+shortoptions[i]+ "\t: " + usage[i];
}
break;
case INT :
{
intvalues[longoptions[i]] = 0;
struct option o = {longoptions[i].c_str(), required_argument, NULL, shortoptions[i]};
longOpts[i] = o;
shortOpts += shortoptions[i];
shortOpts += ":";
optionroles += "\n--"+longoptions[i]+", -"+shortoptions[i]+ " <int>\t: " + usage[i];
}
break;
}
purpose = _purpose;
fullusage = _fullusage;
}
struct option o = {"help", no_argument, NULL, 'h'};
longOpts[noptions] = o;
shortOpts += "h";
if (flagvalues.size() > 0){
synopsis += " [-";
for (map<string, int>::iterator it = flagvalues.begin(); it != flagvalues.end(); it++)
synopsis += shortnames[it->first];
synopsis += "] ";
}
for (map<string, int>::iterator it = intvalues.begin(); it != intvalues.end(); it++){
synopsis += "[-";
synopsis += shortnames[it->first];
synopsis += " <int>] ";
}
for (map<string, double>::iterator it = floatvalues.begin(); it != floatvalues.end(); it++){
synopsis += "[-";
synopsis += shortnames[it->first];
synopsis += " <float>] ";
}
for (map<string, string>::iterator it = stringvalues.begin(); it != stringvalues.end(); it++){
synopsis += "[-";
synopsis += shortnames[it->first];
synopsis += " <string>] ";
}
synopsis += fullusage+"\n";
}
void Args::info(){
for (map<string, int>::iterator it = intvalues.begin(); it != intvalues.end(); it++)
cout << it->first << ": " << it->second <<endl;
for (map<string, double>::iterator it = floatvalues.begin(); it != floatvalues.end(); it++)
cout << it->first << ": " << it->second <<endl;
for (map<string, string>::iterator it = stringvalues.begin(); it != stringvalues.end(); it++)
cout << it->first << ": " << it->second <<endl;
for (map<string, int>::iterator it = flagvalues.begin(); it != flagvalues.end(); it++)
cout << it->first << ": " << it->second <<endl;
cout << nposarg << " positionnal arguments:\n";
for(int i = 0; i<nposarg; i++)
cout << posarg[i] << endl;
}
void Args::printusage(){
cerr << synopsis << endl
<< purpose << endl
<< optionroles << endl;
exit(EXIT_FAILURE);
}
void Args::processarg(int val){
istringstream s;
if (val == '?'){
cerr << "Incorrect option: " << (char)optopt << endl;
printusage();
}
if ( val == 'h'){
cerr << "Usage: " << endl;
printusage();
}
switch (types[val]){
case STRING:
stringvalues[longnames[val]] = optarg;
break;
case FLOAT:
s.str(optarg);
s >> floatvalues[longnames[val]];
break;
case FLAG:
flagvalues[longnames[val]]++;
break;
case INT:
s.str(optarg);
s >> intvalues[longnames[val]];
break;
}
isset[longnames[val]]= true;
}
void Args::processargs(int argc, char **argv){
int opt = 0;
int longIndex = 0;
string name = argv[0];
synopsis = name.substr(name.find_last_of('/')+1) + synopsis;
synopsis = "Synopsis:\n" + synopsis;
opterr = 0;
opt = getopt_long( argc, argv, shortOpts.c_str(), longOpts, &longIndex );
while( opt != -1 ) {
processarg(opt);
opt = getopt_long( argc, argv, shortOpts.c_str(), longOpts, &longIndex );
}
nposarg =argc - optind;
posarg = new string[nposarg];
for (int i = optind; i < argc; i++)
posarg[i - optind] = argv[i];
}
// template <int>
// void Args::get(int prout, string key){
// // return intvalues[key];
// }
//double Args::get (string key){return floatvalues[key];}
//string Args::get (string key){return stringvalues[key];}
//bool Args::get (string key){return flagvalues[key];}
bool Args::operator [](string key){
return isset.find(key)->second;
}
string Args::operator[](int i){
return posarg[i];
}
void Args::def (int val,string key){
map<string,int>::iterator it;
it = intvalues.find(key);
if (it != intvalues.end() && (!isset[key]))
it->second = val;
}
void Args::def (double val, string key){
map<string,double>::iterator it;
it = floatvalues.find(key);
if (it != floatvalues.end() && (!isset[key]))
it->second = val;
}
void Args::def (string val, string key){
map<string,string>::iterator it;
it = stringvalues.find(key);
if (it != stringvalues.end() && (!isset[key]))
it->second = val;
}
// Copyright (C) 2008 Marc Betoule
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.