Commit 2a128c5e authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

+ include

parent e60eea25
/*
* Project: Spherelib
* Date: 09/04
* 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 "CovMat.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>
using namespace std;
/////////////////////////////////////////////////
/// Return a gsl matrix for bin z
/// \param CMat : Input covariance matrices object
/// \param z : bin number
/// \param M : pointor to an allocated gsl matrix
/////////////////////////////////////////////////
template <typename T> void setMatrix (const arr2<T> &CMat, int z, int nbin, gsl_matrix *M)
{
int k = 0;
for (int i=0; i<nbin; i++)
for (int j=0; j<=i; j++)
{
gsl_matrix_set (M, i, j, CMat[k][z]);
gsl_matrix_set (M, j, i, CMat[k++][z]);
}
}
/////////////////////////////////////////////////
/// compute ILC filter :
/// \param CMat : Input covariance matrices object
/// \param Filter : array of filter coefficients
/// \param MixCol : mixing column
/////////////////////////////////////////////////
template <typename T> Status ILC (CovMat<T> CMat, arr2<T> &Filter, vector<double> MixCol)
{
int NbIn = CMat.getNbIn();
int NbBin = CMat.getNbBin();
Filter.alloc(NbIn, NbBin);
gsl_matrix * M, *iM;
gsl_vector *a, *Ma;
gsl_permutation *p;
double aMa;
M = gsl_matrix_alloc (NbIn, NbIn);
iM = gsl_matrix_alloc (NbIn, NbIn);
a = gsl_vector_alloc (NbIn);
Ma = gsl_vector_alloc (NbIn);
p = gsl_permutation_alloc (NbIn);
CBLAS_UPLO_t Uplo=CblasUpper;
for (int i=0; i<NbIn; i++)
gsl_vector_set(a, i, MixCol[i]);
int signum;
const arr2<T> CArr = CMat.getCovMat();
for (long z=0; z<NbBin; z++)
{
setMatrix(CArr, z, NbIn, M);
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);
gsl_blas_ddot (a, Ma, &aMa);
gsl_vector_scale (Ma, (1.0/aMa));
for (int i=0; i<NbIn; i++)
Filter[i][z] = gsl_vector_get (Ma, i);
}
gsl_matrix_free (M);
gsl_matrix_free (iM);
gsl_vector_free (a);
gsl_vector_free (Ma);
gsl_permutation_free (p);
return 0;
}
/////////////////////////////////////////////////
/// compute matched filter :
/// \param CMat : Input covariance matrices object
/// \param Filter : array of filter coefficients
/// \param MixCol : mixing column
/////////////////////////////////////////////////
template <typename T> Status MF (CovMat<T> CMat, arr2<T> &Filter, vector<double> MixCol, arr<T> &noise, arr2<T> &svd)
{
int NbIn = CMat.getNbIn();
int NbBin = CMat.getNbBin();
Filter.alloc(NbIn, NbBin);
gsl_vector * A = gsl_vector_alloc (NbIn);
gsl_vector * AtinvR = gsl_vector_alloc (NbIn);
for (int i=0; i<NbIn; i++)
gsl_vector_set(A, i, MixCol[i]);
gsl_matrix * R = gsl_matrix_alloc (NbIn, NbIn);
gsl_matrix * V = gsl_matrix_alloc (NbIn, NbIn);
gsl_vector * S = gsl_vector_alloc (NbIn);
gsl_vector * work = gsl_vector_alloc (NbIn);
arr2<T> CArr = CMat.getCovMat();
double eps = 2.2204e-16;
for (long z=0; z<NbBin; z++)
{
setMatrix(CArr, z, NbIn, R);
gsl_linalg_SV_decomp (R, V, S, work);
double cond;
int rank = 1;
double result;
cond = gsl_vector_get(S,0)/gsl_vector_get(S,NbIn-1);
for (int i = NbIn - 1; i > 0; i--){
if (gsl_vector_get(S,i) < (gsl_vector_get(S,0)*eps)) gsl_vector_set(S,i,0.);
else{rank = i + 1;break;}
}
gsl_linalg_SV_solve (R, V, S, A, AtinvR);
gsl_blas_ddot (AtinvR , A, &result);
noise[z] = result;
svd[0][z] = cond;
svd[1][z] = rank;
for (int i=0; i<NbIn; i++)
Filter[i][z] = gsl_vector_get (AtinvR, i);
}
gsl_matrix_free (R);
gsl_matrix_free (V);
gsl_vector_free (S);
gsl_vector_free (work);
gsl_vector_free (A);
gsl_vector_free (AtinvR);
return 0;
}
/////////////////////////////////////////////////
/// compute spectral mismatch
/// \param CMat : Input covariance matrices object
/// \param CMat2 : Input covariance matrices object
/// \param wq : number of modes
/////////////////////////////////////////////////
template <typename T>
double kullback (CovMat<T> CMat, CovMat<T> CMat2, arr<T> wq)
{
// init mismatch
double res = 0;
double tmp = 0;
int m = CMat.getNbIn();
int nq = CMat.getNbBin();
// used to store covariance matrices
gsl_matrix *hR, *R, *iR;
hR = gsl_matrix_alloc (m, m);
R = gsl_matrix_alloc (m, m);
iR = gsl_matrix_alloc (m, m);
// used to invert Rmodel
CBLAS_UPLO_t Uplo=CblasLower;
int signum;
gsl_permutation *p;
p = gsl_permutation_alloc (m);
// used to compute matrix-matrix product
CBLAS_SIDE_t Side=CblasLeft;
const arr2<T> CArr = CMat.getCovMat();
const arr2<T> CArr2 = CMat2.getCovMat();
// used to compute svd
gsl_vector_complex *eig;
gsl_eigen_nonsymm_workspace *work;
work = gsl_eigen_nonsymm_alloc (m);
eig = gsl_vector_complex_alloc (m);
// compute mismatch
int q, i, j, k;
int ncross = (int)m*(m+1)/2;
for (q=0; q<nq; q++)
{
tmp = 0;
//init R : 2 cases
// contaminant covariance is not fixed
// cmb contribution already contained in smic->rq
setMatrix(CArr2, q, m, R);
setMatrix(CArr, q, m, hR);
//invert model
gsl_linalg_LU_decomp (R, p, &signum);
gsl_linalg_LU_invert (R, p, iR);
for (i=0; i<m; i++)
for (j=0; j<=i; j++){
gsl_matrix_set(R, i, j, 0) ;
}
// product
// now R contains hR*iR
gsl_blas_dsymm (Side, Uplo, 1, hR, iR, 1, R);
for (i=0; i<m; i++)
for (j=0; j<=i; j++)
gsl_matrix_get(R, i, j);
// svd
gsl_eigen_nonsymm (R, eig, work);
// sum eigen values
double ei ;
for (i=0; i<m; i++){
ei = gsl_matrix_get(R, i, i);
tmp += ei - log(ei);}
tmp -= m;
tmp *= 0.5*wq[q];
res += tmp;
}
gsl_eigen_nonsymm_free (work);
gsl_permutation_free(p);
gsl_vector_complex_free (eig);
gsl_matrix_free (hR);
gsl_matrix_free (R);
gsl_matrix_free (iR);
return res;
}
// Copyright (C) 2007,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
//Purpose:
//Iteratively fills holes in an healpix map by putting in empty pixels
//the mean value of its neighbors
#define MAX(x,y) x>y?x:y
#include <cstdlib>
#include "fitshandle.h"
#include <vector>
#include <cmath>
#include <sstream>
#include "healpix_map.h"
#include "healpix_map_fitsio.h"
using namespace std;
template <typename T> T fillPixNeig(Healpix_Map<T> & m, vector<int> & l){
//vector<int>::iterator icenter;
T maxdiff=0;
int npix = l.size();
vector<T> val(npix);
#pragma omp parallel
{
T lap;
T localmax = 0;
fix_arr< int, 8 >neigpix;
T diff=0;
#pragma omp for
for(int ipix = 0; ipix < npix; ipix++){
lap = 0;
m.neighbors(l[ipix], neigpix);
for(int i = 0; i < 8; i+=2) lap += m[neigpix[i]];
lap -= 4.*m[l[ipix]];
lap /= 5.0;
//if(0 != mean) diff = abs((mean-m[l[ipix]])/mean);
diff = abs(lap);
if (diff>localmax) localmax=diff;
val[ipix]=lap+m[l[ipix]];
}
#pragma omp for
for(int ipix = 0; ipix < npix; ipix++){
m[l[ipix]] = val[ipix];
}
#pragma omp critical
{
if(localmax>maxdiff) maxdiff = localmax;
}
}
return maxdiff;
}
/*
float fillPixDisc(Healpix_Map<float> & m, vector<int> & l, float radius){
vector<int>::iterator icenter, ineig;
vector<int> neigpix;
float mean;
float maxdiff=0;
float diff=0;
for(icenter = l.begin(); icenter < l.end(); icenter++){
mean = 0;
m.query_disc(m.pix2ang(*icenter), radius, neigpix);
for(ineig = neigpix.begin(); ineig < neigpix.end(); ineig++) mean += m[*ineig];
mean /= (float)(neigpix.size());
if(mean != 0) diff = abs((mean-m[*icenter])/mean);
maxdiff = MAX(maxdiff,diff);
m[*icenter] = mean;
}
return maxdiff;
}
*/
//void solvescale(Healpix_Map<float> & mapin, float eps, float radius, int niter, vector<int> & listEPix){
template <typename T> void solvescale(Healpix_Map<T> &mapin, T eps, int niter, vector<int> & listEPix, bool verbose=0){
if (verbose)
cout << "Iter : ";
//Iterative solve
T maxdiff = 2*eps;
for(int iter = 0; (iter != niter) and (maxdiff > eps) ; iter++){
//if(radius <= 0) maxdiff = fillPixNeig(mapin, listEPix);
//else maxdiff = fillPixDisc(mapin, listEPix, radius);
maxdiff = fillPixNeig(mapin, listEPix);
if (verbose)
cout << flush << iter << "(" << maxdiff << ")"<< " . ";
}
if (verbose)
cout << "\n";
}
template <typename T> int fastfiller (Healpix_Map<T> &mapin, Healpix_Map<T> &mask,T eps=1e-6, int niter=-1, int startscale=0, bool verbose=0)
{
//double radius = 0;
//radius = radius/60*degr2rad; //convert radius in radian
//store some generic info on the map
bool swapped = false; // if it was ring ordered first then put it back to ring ordering scheme after computation
T maxvalue = mapin[0]; // use as a scale for the convergence criterion.
vector<int> finestListEPix; // list of pixels to fill
if(startscale == 0) startscale = mapin.Nside();
if(startscale > 256) startscale = 256;
if(startscale < mapin.Nside()){ //save unusefull computation in case of monoscale resolution
//translate the map in Nest ordering scheme
if(mapin.Scheme() == RING){mapin.swap_scheme(); mask.swap_scheme(); swapped=true;}
//Use Healpix_undef rather than 0
for(int i = 0; i < mapin.Npix(); i++){
if(mapin[i] == 0 || approx<T>(mapin[i], Healpix_undef)){
if (mask[i] !=0){
mapin[i] = Healpix_undef;
finestListEPix.push_back(i);
}}
else maxvalue = MAX(maxvalue,mapin[i]);
}
//Compute the coarsest solution
Healpix_Map<T> startmap(startscale, mapin.Scheme(), SET_NSIDE);
startmap.Import_degrade(mapin);
vector<int> listEPix;
for(int i = 0; i < startmap.Npix(); i++){
if(approx<T>(startmap[i], Healpix_undef)){
startmap[i] = 0;
listEPix.push_back(i);
}
}
solvescale(startmap, eps*maxvalue, niter, listEPix, verbose);
//fitshandle outfits;
//outfits.create("!coarstest.fits");
//write_Healpix_map_to_fits (outfits,startmap,FITSUTIL<float>::DTYPE);
//refine this solution
for(int scale = startscale*2; scale < mapin.Nside(); scale*=2){
Healpix_Map<T> mapscale(scale, mapin.Scheme(), SET_NSIDE);
mapscale.Import_degrade(mapin);
//Build the list of empty pixels and initialize them with the coarser solution
listEPix.clear();
for(int i = 0; i < mapscale.Npix(); i++){
if(approx<T>(mapscale[i], Healpix_undef)){
mapscale[i] = startmap[i/4];
listEPix.push_back(i);
}
}
//Solve the finer solution
solvescale(mapscale, eps*maxvalue, niter, listEPix, verbose);
//ostringstream filename;
//filename << "!scale" << scale << ".fits";
//outfits.create(filename.str());
//take this as starting point for the next refinement step
startmap.swap(mapscale);
//write_Healpix_map_to_fits (outfits,startmap,FITSUTIL<float>::DTYPE);
}
//Update the mapin
#pragma omp parallel for
for(int i = 0; i < mapin.Npix(); i++){
if(approx<T>(mapin[i],Healpix_undef)){
mapin[i] = startmap[i/4];
}
}
}
else{
for(int i = 0; i < mapin.Npix(); i++){
if(mapin[i] == 0 || approx<T>(mapin[i], Healpix_undef))
finestListEPix.push_back(i);
else
maxvalue = MAX(maxvalue,mapin[i]);
}
}
//Compute the finest solution
if (verbose)
cout << "eps is " << eps*maxvalue << "\n";
solvescale(mapin, eps*maxvalue, niter, finestListEPix, verbose);
if(swapped) mapin.swap_scheme();
return 0;
}
/*
* Project: Spherelib
* Date: 09/04
* Authors: M. Betoule 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 "MapExt.h"
#include "window.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_linalg.h>
using namespace std;
/////////////////////////////////////////////////
/// Build a mask made of KP zones
/// \param covers : coverage fractions
/// \param Map : output zone map
/// \param KMap : reference intensity map
/////////////////////////////////////////////////
template <typename T> Status KPMask (arr<T> covers, MapExt<T> &Map, const MapExt<T> KMap)
{
int ncovers = covers.size();
if ((ncovers == 1 ) && (covers[0]==1))
{
Map.fill(0);
return 0;
}
// make an histogram of reference map
T thresh;
arr<double> KArr (KMap.Map());
gsl_sort (KArr.begin(), 1, KMap.Npix());
int lmax = 2*KMap.Nside();
AlmExt<T> uKAlm (lmax,lmax);
arr<T> Bl (lmax+1);
GaussBeam (Bl, 120);
PowSpecExt<T> Beam (1, lmax);
Beam.Set (Bl);
MapExt<T> Mask ; Mask.SetNside(KMap.Nside(), RING); Mask.fill(-1);
MapExt<T> dKMap ; dKMap.SetNside(64, RING);
MapExt<T> uKMap ; uKMap.SetNside(KMap.Nside(), RING);
for (int c=0; c<ncovers; c++)
{
thresh = (T) KArr[(int)(ceil(covers[c]*KMap.Npix()))];
dKMap.Import (KMap);
for (int p=0; p<dKMap.Npix(); p++)
if (dKMap[p]<thresh)
dKMap[p]=1;
else
dKMap[p]=0;
uKMap.Import(dKMap);
uKAlm.FromMap(uKMap);
uKAlm.ScaleL (Beam.tt());
uKMap.FromAlm(uKAlm);
for (int p=0; p<uKMap.Npix(); p++)
if ( (uKMap[p]<0.5) && (Mask[p]==-1))
Mask[p] = c;
}
for (int p=0; p<uKMap.Npix(); p++)
if ((Mask[p]==-1))
Mask[p] = ncovers;
Map.Import(Mask);
for (int p=0; p<Map.Npix(); p++)
Map[p] = round(Map[p]);
return 0;
}
/////////////////////////////////////////////////
/// Compute template maps up to order
/// \param Tpl : template map array
/// \param order : regression order (1 or 2)
/////////////////////////////////////////////////
template <typename T> Status getTemplate (MapExt<T> *Map, int order, int nside, Healpix_Ordering_Scheme scheme)
{
(Map[0]).SetNside (nside, scheme);
(Map[0]).fill(1.0);
if (order != 1)
{
(Map[1]).SetNside (nside, scheme);
(Map[2]).SetNside (nside, scheme);
(Map[3]).SetNside (nside, scheme);
for (int pix=0; pix<12*nside*nside; pix++)
{
pointing p = (Map[1]).pix2ang (pix);
(Map[1])[pix] = sin(p.theta)*cos(p.phi);
(Map[2])[pix] = sin(p.theta)*sin(p.phi);
(Map[3])[pix] = cos(p.theta);
}
}
}
/////////////////////////////////////////////////
/// Remove multipole of masked map up to order
/// \param Mask : mask that will be applied
/// \param order : regression order (1 or 2)
/////////////////////////////////////////////////
template <typename T> Status Regress (MapExt<T> &Map, const MapExt<T> &Mask, int order=1, bool verbose=0)
{
// get regression order
int ntpl;
if (order==1)
ntpl = 1;
else {
order= 2;
ntpl = 4;
}
// get template
MapExt<T> Tpl [ntpl];
getTemplate (Tpl, order, Map.Nside() , Map.Scheme());
// mask * template
MapExt<T> MskTpl [ntpl];
for (int i=0; i<ntpl; i++){
MskTpl[i] = Tpl[i];
MskTpl[i] *= Mask;
}
gsl_matrix * PtPM, *iPtPM;
gsl_permutation *p;
int signum;
p = gsl_permutation_alloc (ntpl);
PtPM = gsl_matrix_alloc (ntpl, ntpl);
iPtPM = gsl_matrix_alloc (ntpl, ntpl);
for (int i=0; i<ntpl; i++){
for (int j=i; j<ntpl; j++)
{
MapExt<T> tmpMap = Tpl[i];
tmpMap *= MskTpl[j]; // template2 * mask
T val = tmpMap.sum();
gsl_matrix_set (PtPM, i, j, val);
gsl_matrix_set (PtPM, j, i, val);
}
}
if (verbose){
cout << "PtPM is : \n";
for (int i=0; i<ntpl; i++){
for (int j=i; j<ntpl; j++)
{
cout << " " << gsl_matrix_get (PtPM, i, j);
cout << " " << gsl_matrix_get (PtPM, j, i);