Commit 9a3dc517 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

retrieve attempt

parent 010e916d
......@@ -74,7 +74,9 @@ template <typename T> Status ILC (CovMat<T> CMat, arr2<T> &Filter, vector<double
for (long z=0; z<NbBin; z++)
{
setMatrix(CArr, z, NbIn, M);
gsl_matrix_set_zero (iM);
gsl_vector_set_zero (Ma);
gsl_linalg_LU_decomp (M, p, &signum);
gsl_linalg_LU_invert (M, p, iM);
gsl_blas_dsymv (Uplo, 1.0, iM, a, 1.0, Ma);
......
......@@ -58,6 +58,20 @@ void print_time (time_t start, time_t startu, time_t end, time_t endu){
double elapsed_useconds = endu-startu;
fprintf(stdout, "elapsed time %lf s\n", elapsed_seconds + elapsed_useconds/1000000.0);
}
void rm_tmp_file (string str){
char tplt[100];
sprintf(tplt, "rm -f %s", str.c_str());
system(tplt);
}
void rm_lst_tmp (int nfile, string** arr){
char tplt[100];
for (int i=0; i<nfile; i++) {
sprintf(tplt, "rm -r %s", (*arr)[i].c_str());
system(tplt);
}
}
string get_tmp_file (string str){
string tempdir = getenv("TMPDIR");
......@@ -1092,6 +1106,97 @@ template <typename T> int map2alm (string infile, string outfile, int lmax, int
}
/////////////////////////////////////////////////
/// Compute spherical harmonic coefficients from Healpix map
/// \param infile : map file name
/// \param outfile : out alm file name
/// \param lmax : maximum order of analysis
/// \param nside : needlet resolution
/// \param niter : iteration number
/// \param filterfile : needlet window function
template <typename T> int map2need (string infile, string outfile, int lmax, int nside, int niter,string filterfile)
{
fitshandle out;
AlmExt <T> alm;
MapExt <T> map;
read_Healpix_map_from_fits(infile, map, 1, 2);
alm.FromMap (map, lmax, lmax, niter); // get alm
out.create (outfile);// save them to file
// init filter reading specfile
arr<T> filter (lmax);
readfitscol(filterfile, filter, 1);
alm.ScaleL (filter);
MapExt <T> need;
need.FromAlm (alm, nside);// compute Healpix map
write_Healpix_map_to_fits (out, need,planckType<T>());// save it to file
return 0;
}
/////////////////////////////////////////////////
/// Compute spherical harmonic coefficients from Healpix map
/// \param infile : map file name
/// \param outfile : out alm file name
/// \param lmax : maximum order of analysis
/// \param niter : iteration number
/// \param filterfile : needlet window function
template <typename T> int need2alm (string infile, string outfile, int lmax, int niter,string filterfile)
{
fitshandle out;
AlmExt <T> alm;
MapExt <T> map;
read_Healpix_map_from_fits(infile, map, 1, 2);
alm.FromMap (map, lmax, lmax, niter); // get alm
out.create (outfile);// save them to file
// init filter reading specfile
arr<T> filter (lmax);
readfitscol(filterfile, filter, 1);
alm.ScaleL (filter);
write_Alm_to_fits (out,alm,alm.Lmax(),alm.Mmax(),planckType<T>());
return 0;
}
/////////////////////////////////////////////////
/// Compute spherical harmonic coefficients from Healpix map
/// \param infile : alm file name
/// \param outfile : out needlet file name
/// \param lmax : maximum order of analysis
/// \param nside : needlet resolution
/// \param filterfile : needlet window function
template <typename T> int alm2need (string infile, string outfile, int lmax, int nside, string filterfile)
{
fitshandle out;
AlmExt <T> alm;
MapExt <T> map;
read_Alm_from_fits (infile, alm, lmax, lmax);
out.create (outfile);// save them to file
// init filter reading specfile
arr<T> filter (lmax);
readfitscol(filterfile, filter, 1);
alm.ScaleL (filter);
MapExt <T> need;
need.FromAlm (alm, nside);// compute Healpix map
write_Healpix_map_to_fits (out, need,planckType<T>());// save it to file
return 0;
}
/////////////////////////////////////////////////
/// Generate mask
/// \param covers : coverage fractions
......
......@@ -131,7 +131,7 @@ class CovMat
// get number of zones
T2 min , max;
Mask.minmax(min,max);
m_NbBin = (long) max+1;
m_NbBin = (long) max;
// set size of covmat
m_NbCross = m_NbIn*(m_NbIn+1)/2;
......@@ -165,8 +165,8 @@ class CovMat
#pragma omp for
for (long pix=0; pix<Mask.Npix() ; pix++)
{
tempval[(int)Mask[pix]] += (Maps[m1])[pix] * (Maps[m2])[pix];// sum product
tempwei[(int)Mask[pix]] += 1;
tempval[(int)Mask[pix]-1] += (Maps[m1])[pix] * (Maps[m2])[pix];// sum product
tempwei[(int)Mask[pix]-1] += 1;
}
#pragma omp critical
{
......@@ -206,7 +206,7 @@ class CovMat
// get number of zones
T2 min , max;
Mask.minmax(min,max);
m_NbBin = (long) max+1;
m_NbBin = (long) max;
// set size of covmat
m_NbCross = m_NbIn*(m_NbIn+1)/2;
......@@ -235,8 +235,8 @@ class CovMat
#pragma omp for
for (long pix=0; pix<Mask.Npix() ; pix++)
{
tempval[(int)Mask[pix]] += (Maps[m1])[pix];
tempwei[(int)Mask[pix]] += 1;
tempval[(int)Mask[pix]-1] += (Maps[m1])[pix];
tempwei[(int)Mask[pix]-1] += 1;
}
#pragma omp critical
{
......@@ -275,8 +275,8 @@ class CovMat
#pragma omp for
for (long pix=0; pix<Mask.Npix() ; pix++)
{
tempval[(int)Mask[pix]] += ((Maps[m1])[pix]-Vmean[m1][(int)Mask[pix]]) * ((Maps[m2])[pix]-Vmean[m2][(int)Mask[pix]]);// sum product
tempwei[(int)Mask[pix]] += 1;
tempval[(int)Mask[pix]-1] += ((Maps[m1])[pix]-Vmean[m1][(int)Mask[pix]-1]) * ((Maps[m2])[pix]-Vmean[m2][(int)Mask[pix]-1]);// sum product
tempwei[(int)Mask[pix]-1] += 1;
}
#pragma omp critical
{
......
......@@ -131,10 +131,10 @@ class MapExt : public Healpix_Map<T>
CD_iCheck( Conformable(Mask) );
T2 min, max;
Mask.minmax(min,max);
long nzone = (long) max+1;
long nzone = (long) max;
// for each pixel if zone z
for (int pix=0; pix<Mask.Npix(); pix++)
(*this)[pix] += beta[(int)Mask[pix]]*MapB[pix];
(*this)[pix] += beta[(int)Mask[pix]-1]*MapB[pix];
CD_returnOK;
......@@ -168,7 +168,7 @@ class MapExt : public Healpix_Map<T>
CD_iCheck( Conformable(Mask) );
T2 min, max;
Mask.minmax(min,max);
long nzone = (long) max+1;
long nzone = (long) max;
Mean.alloc(nzone);
Mean.fill(0);
arr<int>npix(nzone);
......@@ -179,8 +179,8 @@ class MapExt : public Healpix_Map<T>
#pragma omp for
// for each pixel if zone z
for (int pix=0; pix<Mask.Npix(); pix++){
Mean[Mask[pix]] += (*this)[pix];
npix[Mask[pix]]++;}
Mean[Mask[pix]-1] += (*this)[pix];
npix[Mask[pix]-1]++;}
}
for (int z=0; z<nzone; z++)
Mean[z] = Mean[z]/npix[z];
......@@ -203,11 +203,11 @@ class MapExt : public Healpix_Map<T>
CD_iCheck( Conformable(Mask) );
T2 min, max;
Mask.minmax(min,max);
long nzone = (long) max+1;
long nzone = (long) max;
// for each pixel if zone z
for (int pix=0; pix<Mask.Npix(); pix++)
for(long Idx=0;Idx < NbMap; Idx++)
(*this)[pix] += ListCoef[Idx][Mask[pix]]* (*ListMap[Idx])[pix];
(*this)[pix] += ListCoef[Idx][Mask[pix]-1]* (*ListMap[Idx])[pix];
CD_returnOK;
}
......
......@@ -18,7 +18,16 @@ 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,"."] )
AddOption('--prefix',
dest='prefix',
type='string',
nargs=1,
action='store',
metavar='DIR',
default='/usr/local',
help='installation prefix')
opt = Environment(PREFIX = GetOption('prefix'), 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)')
......@@ -36,5 +45,5 @@ 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')
opt.Install('$PREFIX/bin',opt.Program(cfiles[p]) )
opt.Alias('install', '$PREFIX/bin')
......@@ -54,11 +54,11 @@ int main(int argc ,const char** argv)
// 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 ");
string longopt[] = {"almopt", "n_iter", "nside", "pixel_based", "zonefile"};
ArgType type[] = {FLAG, INT, INT, FLAG, STRING};
string usage[] = {"Inputs are alm FITS file.", "Number of iteration of anafast.", "Resolution of the output map.", "Pixel based analysis", "Pixel zone definition"};
char shortopt[] = "ainpz";
Args args = Args(5, longopt, shortopt, type, usage,"Perform needlet ILC.", "outfile l0 J zonefile1 lmax1 ... zonefileJ lmaxJ N infile1 coeff1 ... infileN coeffN \n -p outfile -z zonefile N infile1 coeff1 ... infileN coeffN ");
// parameter file
......@@ -110,28 +110,47 @@ int main(int argc ,const char** argv)
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++]);}
args.get(n_iter,"n_iter");
args.get(nside,"nside");
pixel_based = args["pixel_based"];
if (pixel_based){
string zfile = "";
args.get(zfile,"zonefile");
if (zfile.size()>1){
zonefile = new string[1];
zonefile[0] = zfile;
nbands = 1;}
else
nbands = 0;
int idx=0;
nbin = atoi(args[idx++].c_str());
infile = new string[nbin];
mixcol.alloc(nbin);
for (int i=0; i<nbin; i++){
infile[i] = args[idx++];
mixcol[i] = atof(args[idx++].c_str());}
}
else{
int idx = 1;
outfile = args[idx++];
int l0 = atoi(args[idx++].c_str());
nbands = atoi(args[idx++].c_str());
bands.alloc(nbands+2);
zonefile = new string[nbands];
bands[0] = 0;
bands[1] = l0;
for (int i=0; i<nbands; i++){
zonefile[i] = args[idx++];
bands[i+2] = atoi(args[idx++].c_str());}
nbin = atoi(args[idx++].c_str());
infile = new string[nbin];
mixcol.alloc(nbin);
for (int i=0; i<nbin; i++){
infile[i] = args[idx++];
mixcol[i] = atof(args[idx++].c_str());}
}
}
// preprocessing
fitshandle inp;
if (!almopt){
......@@ -166,31 +185,22 @@ int main(int argc ,const char** argv)
if (pixel_based){
string statfile = get_tmp_file("stat");
string ilcfile = get_tmp_file("ilc");
fprintf(stdout, "Compute second order stats: ");
gettimeofday(&start, NULL);
fprintf(stdout, "Compute second order stats\n");
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);
infile = mapfile;}
map2stat<double> (infile, nbin, statfile, zonefile[0]);
fprintf(stdout, "Compute ILC filter\n");
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);
fprintf(stdout, "Apply ILC filter\n ");
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);
rm_tmp_file(statfile);
rm_tmp_file(ilcfile);
}
else{
string *winfile;
get_lst_tmp (nbands, "window", &winfile);
......@@ -202,45 +212,33 @@ int main(int argc ,const char** argv)
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);
fprintf(stdout, "Perform analysis of map %d and scale %d. Lmax is %d, Nside is %d\n ", m, j, bands[j+2], nsides[j]);
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);
alm2need<double> (infile[m], mapfile[m], bands[j+2], nsides[j], winfile[j]);
else
map2need<double> (infile[m], mapfile[m], bands[j+2], nsides[j], n_iter,winfile[j]);
}
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);
fprintf(stdout, "Compute ILC filter for scale %d\n ", j);
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);
fprintf(stdout, "Apply ILC filter for scale %d\n ", j);
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);
need2alm<double> (cmapfile, scalefile[j], bands[j+2], n_iter, winfile[j]);
rm_tmp_file(statfile);
rm_tmp_file(ilcfile);
rm_tmp_file(cmapfile);
rm_lst_tmp(nbin, &mapfile);
}
string almfile = get_tmp_file("alm");
fprintf(stdout, "Combine scales: ");
gettimeofday(&start, NULL);
fprintf(stdout, "Combine scales\n ");
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);
rm_lst_tmp(nbands ,&scalefile);
rm_lst_tmp(nbands, &winfile);
rm_tmp_file(almfile);
}
return 0;
......
......@@ -85,7 +85,7 @@ def map2cov (maptab, zone, reg=0):
nin = shape(maptab)[0]
npix = max(zone)+1
npix = max(zone)
ncross = nin*(nin+1)/2
covtab = zeros((ncross, npix))
nmode = zeros((npix))
......
......@@ -87,10 +87,10 @@ def build(bld):
cxx_atlas = bld(
features = ['cxx', 'cshlib'],
source = ['../lib/src/catalog.cpp', '../lib/src/profile.cpp', '../lib/src/tangentPlane.cpp', '../lib/src/maptools.cpp', '../lib/src/atlasmanager.cpp'],
source = ['../lib/atlas/catalog.cpp', '../lib/atlas/profile.cpp', '../lib/atlas/tangentPlane.cpp', '../lib/atlas/maptools.cpp', '../lib/atlas/atlasmanager.cpp'],
target = 'libatlas',
vnum = get_version(),
includes = ['../lib/src', '../healpix/'+target+'/include'],
includes = ['../lib/atlas', '../healpix/'+target+'/include'],
defines = [hp_data],
ccflags = bld.env['CXXFLAGS'],
lib = libs+['fftw3', 'm'],
......@@ -103,7 +103,7 @@ def build(bld):
swig_spherelib_map = bld(
features = 'cxx cshlib pyext',
includes = ['../lib/src', '../healpix/'+target+'/include','spherelib',
includes = ['../lib/atlas','../lib/src', '../healpix/'+target+'/include','spherelib',
'../include', '../lib/src'],
source = ['spherelib/spherelib_map.i','spherelib/spherelib_map.cpp'],
target = '_spherelib_map',
......@@ -155,5 +155,5 @@ def build(bld):
obj.install_path = '${PYTHONDIR}/spherelib'
## TODO !!
print('→ building icosahedron')
os.system("gfortran ../include/icosahedron.f -o "+'${PREFIX}'+"/bin/icosahedron")
#print('→ building icosahedron')
#os.system("gfortran ../include/icosahedron.f -o "+'${PREFIX}'+"/bin/icosahedron")
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