Skip to content

Commit 1855928

Browse files
damongedamonge
committed
implemented, not tested
1 parent 5df0743 commit 1855928

File tree

6 files changed

+295
-0
lines changed

6 files changed

+295
-0
lines changed

‎pymaster/namaster.i‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -489,6 +489,16 @@ nmt_workspace *comp_coupling_matrix(int spin1,int spin2,
489489
norm_type, w2);
490490
}
491491

492+
void comp_general_coupling_matrix(int s1, int s2, int n1, int n2, int lmax,
493+
int nell4,double *f_ell,
494+
double *dout,int ndout)
495+
{
496+
asserting(nell4==lmax+1);
497+
asserting(ndout==nell4*nell4);
498+
memset(dout,0,ndout*sizeof(double));
499+
nmt_compute_general_coupling_matrix(lmax,f_ell,s1,s2,n1,n2,dout);
500+
}
501+
492502
nmt_workspace_flat *comp_coupling_matrix_flat(nmt_field_flat *fl1,nmt_field_flat *fl2,
493503
nmt_binning_scheme_flat *bin,
494504
double lmn_x,double lmx_x,double lmn_y,double lmx_y,

‎pymaster/namaster_wrap.c‎

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3456,6 +3456,16 @@ nmt_workspace *comp_coupling_matrix(int spin1,int spin2,
34563456
norm_type, w2);
34573457
}
34583458

3459+
void comp_general_coupling_matrix(int s1, int s2, int n1, int n2, int lmax,
3460+
int nell4,double *f_ell,
3461+
double *dout,int ndout)
3462+
{
3463+
asserting(nell4==lmax+1);
3464+
asserting(ndout==nell4*nell4);
3465+
memset(dout,0,ndout*sizeof(double));
3466+
nmt_compute_general_coupling_matrix(lmax,f_ell,s1,s2,n1,n2,dout);
3467+
}
3468+
34593469
nmt_workspace_flat *comp_coupling_matrix_flat(nmt_field_flat *fl1,nmt_field_flat *fl2,
34603470
nmt_binning_scheme_flat *bin,
34613471
double lmn_x,double lmx_x,double lmn_y,double lmx_y,
@@ -11628,6 +11638,75 @@ SWIGINTERN PyObject *_wrap_compute_coupling_matrix_anisotropic(PyObject *SWIGUNU
1162811638
}
1162911639

1163011640

11641+
SWIGINTERN PyObject *_wrap_compute_general_coupling_matrix(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
11642+
PyObject *resultobj = 0;
11643+
int arg1 ;
11644+
flouble *arg2 = (flouble *) 0 ;
11645+
int arg3 ;
11646+
int arg4 ;
11647+
int arg5 ;
11648+
int arg6 ;
11649+
flouble *arg7 = (flouble *) 0 ;
11650+
int val1 ;
11651+
int ecode1 = 0 ;
11652+
void *argp2 = 0 ;
11653+
int res2 = 0 ;
11654+
int val3 ;
11655+
int ecode3 = 0 ;
11656+
int val4 ;
11657+
int ecode4 = 0 ;
11658+
int val5 ;
11659+
int ecode5 = 0 ;
11660+
int val6 ;
11661+
int ecode6 = 0 ;
11662+
void *argp7 = 0 ;
11663+
int res7 = 0 ;
11664+
PyObject *swig_obj[7] ;
11665+
11666+
if (!SWIG_Python_UnpackTuple(args, "compute_general_coupling_matrix", 7, 7, swig_obj)) SWIG_fail;
11667+
ecode1 = SWIG_AsVal_int(swig_obj[0], &val1);
11668+
if (!SWIG_IsOK(ecode1)) {
11669+
SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "compute_general_coupling_matrix" "', argument " "1"" of type '" "int""'");
11670+
}
11671+
arg1 = (int)(val1);
11672+
res2 = SWIG_ConvertPtr(swig_obj[1], &argp2,SWIGTYPE_p_double, 0 | 0 );
11673+
if (!SWIG_IsOK(res2)) {
11674+
SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "compute_general_coupling_matrix" "', argument " "2"" of type '" "flouble *""'");
11675+
}
11676+
arg2 = (flouble *)(argp2);
11677+
ecode3 = SWIG_AsVal_int(swig_obj[2], &val3);
11678+
if (!SWIG_IsOK(ecode3)) {
11679+
SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "compute_general_coupling_matrix" "', argument " "3"" of type '" "int""'");
11680+
}
11681+
arg3 = (int)(val3);
11682+
ecode4 = SWIG_AsVal_int(swig_obj[3], &val4);
11683+
if (!SWIG_IsOK(ecode4)) {
11684+
SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "compute_general_coupling_matrix" "', argument " "4"" of type '" "int""'");
11685+
}
11686+
arg4 = (int)(val4);
11687+
ecode5 = SWIG_AsVal_int(swig_obj[4], &val5);
11688+
if (!SWIG_IsOK(ecode5)) {
11689+
SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "compute_general_coupling_matrix" "', argument " "5"" of type '" "int""'");
11690+
}
11691+
arg5 = (int)(val5);
11692+
ecode6 = SWIG_AsVal_int(swig_obj[5], &val6);
11693+
if (!SWIG_IsOK(ecode6)) {
11694+
SWIG_exception_fail(SWIG_ArgError(ecode6), "in method '" "compute_general_coupling_matrix" "', argument " "6"" of type '" "int""'");
11695+
}
11696+
arg6 = (int)(val6);
11697+
res7 = SWIG_ConvertPtr(swig_obj[6], &argp7,SWIGTYPE_p_double, 0 | 0 );
11698+
if (!SWIG_IsOK(res7)) {
11699+
SWIG_exception_fail(SWIG_ArgError(res7), "in method '" "compute_general_coupling_matrix" "', argument " "7"" of type '" "flouble *""'");
11700+
}
11701+
arg7 = (flouble *)(argp7);
11702+
nmt_compute_general_coupling_matrix(arg1,arg2,arg3,arg4,arg5,arg6,arg7);
11703+
resultobj = SWIG_Py_Void();
11704+
return resultobj;
11705+
fail:
11706+
return NULL;
11707+
}
11708+
11709+
1163111710
SWIGINTERN PyObject *_wrap_compute_coupling_matrix(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
1163211711
PyObject *resultobj = 0;
1163311712
int arg1 ;
@@ -16756,6 +16835,116 @@ SWIGINTERN PyObject *_wrap_comp_coupling_matrix(PyObject *SWIGUNUSEDPARM(self),
1675616835
}
1675716836

1675816837

16838+
SWIGINTERN PyObject *_wrap_comp_general_coupling_matrix(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
16839+
PyObject *resultobj = 0;
16840+
int arg1 ;
16841+
int arg2 ;
16842+
int arg3 ;
16843+
int arg4 ;
16844+
int arg5 ;
16845+
int arg6 ;
16846+
double *arg7 = (double *) 0 ;
16847+
double *arg8 = (double *) 0 ;
16848+
int arg9 ;
16849+
int val1 ;
16850+
int ecode1 = 0 ;
16851+
int val2 ;
16852+
int ecode2 = 0 ;
16853+
int val3 ;
16854+
int ecode3 = 0 ;
16855+
int val4 ;
16856+
int ecode4 = 0 ;
16857+
int val5 ;
16858+
int ecode5 = 0 ;
16859+
PyArrayObject *array6 = NULL ;
16860+
int is_new_object6 = 0 ;
16861+
PyObject *array8 = NULL ;
16862+
PyObject *swig_obj[7] ;
16863+
16864+
if (!SWIG_Python_UnpackTuple(args, "comp_general_coupling_matrix", 7, 7, swig_obj)) SWIG_fail;
16865+
ecode1 = SWIG_AsVal_int(swig_obj[0], &val1);
16866+
if (!SWIG_IsOK(ecode1)) {
16867+
SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "comp_general_coupling_matrix" "', argument " "1"" of type '" "int""'");
16868+
}
16869+
arg1 = (int)(val1);
16870+
ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
16871+
if (!SWIG_IsOK(ecode2)) {
16872+
SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "comp_general_coupling_matrix" "', argument " "2"" of type '" "int""'");
16873+
}
16874+
arg2 = (int)(val2);
16875+
ecode3 = SWIG_AsVal_int(swig_obj[2], &val3);
16876+
if (!SWIG_IsOK(ecode3)) {
16877+
SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "comp_general_coupling_matrix" "', argument " "3"" of type '" "int""'");
16878+
}
16879+
arg3 = (int)(val3);
16880+
ecode4 = SWIG_AsVal_int(swig_obj[3], &val4);
16881+
if (!SWIG_IsOK(ecode4)) {
16882+
SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "comp_general_coupling_matrix" "', argument " "4"" of type '" "int""'");
16883+
}
16884+
arg4 = (int)(val4);
16885+
ecode5 = SWIG_AsVal_int(swig_obj[4], &val5);
16886+
if (!SWIG_IsOK(ecode5)) {
16887+
SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "comp_general_coupling_matrix" "', argument " "5"" of type '" "int""'");
16888+
}
16889+
arg5 = (int)(val5);
16890+
{
16891+
npy_intp size[1] = {
16892+
-1
16893+
};
16894+
array6 = obj_to_array_contiguous_allow_conversion(swig_obj[5],
16895+
NPY_DOUBLE,
16896+
&is_new_object6);
16897+
if (!array6 || !require_dimensions(array6, 1) ||
16898+
!require_size(array6, size, 1)) SWIG_fail;
16899+
arg6 = (int) array_size(array6,0);
16900+
arg7 = (double*) array_data(array6);
16901+
}
16902+
{
16903+
npy_intp dims[1];
16904+
if (!PyInt_Check(swig_obj[6]))
16905+
{
16906+
const char* typestring = pytype_string(swig_obj[6]);
16907+
PyErr_Format(PyExc_TypeError,
16908+
"Int dimension expected. '%s' given.",
16909+
typestring);
16910+
SWIG_fail;
16911+
}
16912+
arg9 = (int) PyInt_AsLong(swig_obj[6]);
16913+
dims[0] = (npy_intp) arg9;
16914+
array8 = PyArray_SimpleNew(1, dims, NPY_DOUBLE);
16915+
if (!array8) SWIG_fail;
16916+
arg8 = (double*) array_data(array8);
16917+
}
16918+
{
16919+
try {
16920+
comp_general_coupling_matrix(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9);
16921+
}
16922+
finally {
16923+
SWIG_exception(SWIG_RuntimeError,nmt_error_message);
16924+
}
16925+
}
16926+
resultobj = SWIG_Py_Void();
16927+
{
16928+
resultobj = SWIG_Python_AppendOutput(resultobj,(PyObject*)array8);
16929+
}
16930+
{
16931+
if (is_new_object6 && array6)
16932+
{
16933+
Py_DECREF(array6);
16934+
}
16935+
}
16936+
return resultobj;
16937+
fail:
16938+
{
16939+
if (is_new_object6 && array6)
16940+
{
16941+
Py_DECREF(array6);
16942+
}
16943+
}
16944+
return NULL;
16945+
}
16946+
16947+
1675916948
SWIGINTERN PyObject *_wrap_comp_coupling_matrix_flat(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
1676016949
PyObject *resultobj = 0;
1676116950
nmt_field_flat *arg1 = (nmt_field_flat *) 0 ;
@@ -19356,6 +19545,7 @@ static PyMethodDef SwigMethods[] = {
1935619545
{ "compute_master_coefficients", _wrap_compute_master_coefficients, METH_VARARGS, NULL},
1935719546
{ "master_calculator_free", _wrap_master_calculator_free, METH_O, NULL},
1935819547
{ "compute_coupling_matrix_anisotropic", _wrap_compute_coupling_matrix_anisotropic, METH_VARARGS, NULL},
19548+
{ "compute_general_coupling_matrix", _wrap_compute_general_coupling_matrix, METH_VARARGS, NULL},
1935919549
{ "compute_coupling_matrix", _wrap_compute_coupling_matrix, METH_VARARGS, NULL},
1936019550
{ "update_coupling_matrix", _wrap_update_coupling_matrix, METH_VARARGS, NULL},
1936119551
{ "workspace_update_binning", _wrap_workspace_update_binning, METH_VARARGS, NULL},
@@ -19453,6 +19643,7 @@ static PyMethodDef SwigMethods[] = {
1945319643
{ "synfast_new_flat", _wrap_synfast_new_flat, METH_VARARGS, NULL},
1945419644
{ "comp_coupling_matrix_anisotropic", _wrap_comp_coupling_matrix_anisotropic, METH_VARARGS, NULL},
1945519645
{ "comp_coupling_matrix", _wrap_comp_coupling_matrix, METH_VARARGS, NULL},
19646+
{ "comp_general_coupling_matrix", _wrap_comp_general_coupling_matrix, METH_VARARGS, NULL},
1945619647
{ "comp_coupling_matrix_flat", _wrap_comp_coupling_matrix_flat, METH_VARARGS, NULL},
1945719648
{ "read_workspace", _wrap_read_workspace, METH_VARARGS, NULL},
1945819649
{ "write_workspace", _wrap_write_workspace, METH_VARARGS, NULL},

‎pymaster/nmtlib.py‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -346,6 +346,9 @@ def master_calculator_free(c):
346346
def compute_coupling_matrix_anisotropic(spin1, spin2, mask_aniso_1, mask_aniso_2, lmax, lmax_mask, pcl_masks_00, pcl_masks_0e, pcl_masks_e0, pcl_masks_0b, pcl_masks_b0, pcl_masks_ee, pcl_masks_eb, pcl_masks_be, pcl_masks_bb, beam1, beam2, bin, norm_type, w2):
347347
return _nmtlib.compute_coupling_matrix_anisotropic(spin1, spin2, mask_aniso_1, mask_aniso_2, lmax, lmax_mask, pcl_masks_00, pcl_masks_0e, pcl_masks_e0, pcl_masks_0b, pcl_masks_b0, pcl_masks_ee, pcl_masks_eb, pcl_masks_be, pcl_masks_bb, beam1, beam2, bin, norm_type, w2)
348348

349+
def compute_general_coupling_matrix(lmax, pcl_mask, s1, s2, n1, n2, xi_out):
350+
return _nmtlib.compute_general_coupling_matrix(lmax, pcl_mask, s1, s2, n1, n2, xi_out)
351+
349352
def compute_coupling_matrix(spin1, spin2, lmax, lmax_mask, pure_e1, pure_b1, pure_e2, pure_b2, pcl_masks, beam1, beam2, bin, is_teb, l_toeplitz, l_exact, dl_band, norm_type, w2):
350353
return _nmtlib.compute_coupling_matrix(spin1, spin2, lmax, lmax_mask, pure_e1, pure_b1, pure_e2, pure_b2, pcl_masks, beam1, beam2, bin, is_teb, l_toeplitz, l_exact, dl_band, norm_type, w2)
351354

@@ -535,6 +538,9 @@ def comp_coupling_matrix_anisotropic(spin1, spin2, aniso1, aniso2, lmax, lmax_ma
535538
def comp_coupling_matrix(spin1, spin2, lmax, lmax_mask, pure_e_1, pure_b_1, pure_e_2, pure_b_2, norm_type, w2, nlb1, nlb2, nell4, bin, is_teb, l_toeplitz, l_exact, dl_band):
536539
return _nmtlib.comp_coupling_matrix(spin1, spin2, lmax, lmax_mask, pure_e_1, pure_b_1, pure_e_2, pure_b_2, norm_type, w2, nlb1, nlb2, nell4, bin, is_teb, l_toeplitz, l_exact, dl_band)
537540

541+
def comp_general_coupling_matrix(s1, s2, n1, n2, lmax, nell4, dout):
542+
return _nmtlib.comp_general_coupling_matrix(s1, s2, n1, n2, lmax, nell4, dout)
543+
538544
def comp_coupling_matrix_flat(fl1, fl2, bin, lmn_x, lmx_x, lmn_y, lmx_y, is_teb):
539545
return _nmtlib.comp_coupling_matrix_flat(fl1, fl2, bin, lmn_x, lmx_x, lmn_y, lmx_y, is_teb)
540546

‎pymaster/workspaces.py‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1169,3 +1169,17 @@ def compute_full_master_flat(f1, f2, b, cl_noise=None, cl_guess=None,
11691169
clout = np.reshape(cl1d, [len(cln), b.bin.n_bands])
11701170

11711171
return clout
1172+
1173+
1174+
def get_general_coupling_matrix(pcl_mask, s1, s2, n1, n2):
1175+
""" Returns a general mode-coupling matrix
1176+
TODO
1177+
"""
1178+
1179+
lmax = len(pcl_mask)-1
1180+
xi = lib.comp_general_coupling_matrix(
1181+
int(s1), int(s2), int(n1),
1182+
int(n2), int(lmax),
1183+
pcl_mask, int((lmax+1)**2))
1184+
xi = xi.reshape([lmax+1, lmax+1])
1185+
return xi

‎src/namaster.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -756,6 +756,12 @@ nmt_workspace *nmt_compute_coupling_matrix_anisotropic(int spin1, int spin2,
756756
nmt_binning_scheme *bin,
757757
int norm_type,flouble w2);
758758

759+
void nmt_compute_general_coupling_matrix(int lmax,
760+
flouble *pcl_mask,
761+
int s1, int s2,
762+
int n1, int n2,
763+
flouble *xi_out);
764+
759765
nmt_workspace *nmt_compute_coupling_matrix(int spin1,int spin2,
760766
int lmax, int lmax_mask,
761767
int pure_e1,int pure_b1,

‎src/nmt_master.c‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -994,6 +994,74 @@ nmt_master_calculator *nmt_compute_master_coefficients(int lmax, int lmax_mask,
994994
return c;
995995
}
996996

997+
void nmt_compute_general_coupling_matrix(int lmax,
998+
flouble *pcl_mask,
999+
int s1, int s2,
1000+
int n1, int n2,
1001+
flouble *xi_out)
1002+
{
1003+
int nls=lmax+1;
1004+
int sign = ((n1+n2) & 1) ? -1 : 1;
1005+
1006+
#pragma omp parallel default(none) \
1007+
shared(lmax, pcl_mask, s1, s2, n1, n2, xi_out, nls, sign)
1008+
{
1009+
int same_sn=(s1 == s2) && (n1 == n2);
1010+
int ll2,ll3,icc;
1011+
int lstart=NMT_MAX(s1, s2);
1012+
flouble *wl_mask=my_malloc((lmax+1)*sizeof(flouble));
1013+
double *wigner_sn1=NULL,*wigner_sn2=NULL;
1014+
wigner_sn1=my_malloc(2*(lmax+1)*sizeof(double));
1015+
if(same_sn)
1016+
wigner_sn2=wigner_sn1;
1017+
else
1018+
wigner_sn2=my_malloc(2*(lmax+1)*sizeof(double));
1019+
1020+
for(ll2=0;ll2<=lmax;ll2++)
1021+
wl_mask[ll2]=pcl_mask[ll2]*(2*ll2+1)/(4*M_PI);
1022+
1023+
#pragma omp for schedule(dynamic)
1024+
for(ll2=lstart;ll2<=lmax;ll2++) {
1025+
for(ll3=lstart;ll3<=lmax;ll3++) {
1026+
int l1,lmin_here,lmax_here;
1027+
int lmin_sn1=0,lmax_sn1=2*(lmax+1)+1;
1028+
int lmin_sn2=0,lmax_sn2=2*(lmax+1)+1;
1029+
int index=ll3+(lmax+1)*ll2;
1030+
lmin_here=abs(ll2-ll3);
1031+
lmax_here=ll2+ll3;
1032+
1033+
1034+
drc3jj(ll2,ll3,n1,-s1,&lmin_sn1,&lmax_sn1,wigner_sn1,2*(lmax+1));
1035+
if(same_sn) {
1036+
wigner_sn2=wigner_sn1;
1037+
lmin_sn2=lmin_sn1;
1038+
lmax_sn2=lmax_sn1;
1039+
}
1040+
else
1041+
drc3jj(ll2,ll3,n2,-s2,&lmin_sn2,&lmax_sn2,wigner_sn2,2*(lmax+1));
1042+
1043+
for(l1=lmin_here;l1<=lmax_here;l1++) {
1044+
if(l1<=lmax) {
1045+
int jsn1=l1-lmin_sn1;
1046+
int jsn2=l1-lmin_sn2;
1047+
flouble wsn1=0,wsn2=0;
1048+
wsn1=jsn1 < 0 ? 0 : wigner_sn1[jsn1];
1049+
wsn2=jsn2 < 0 ? 0 : wigner_sn2[jsn2];
1050+
//if(!((l1+ll2+ll3) & 1)) //Even sum
1051+
//if((l1+ll2+ll3) & 1) //Even sum
1052+
xi_out[index] += wl_mask[l1]*wsn1*wsn2;
1053+
}
1054+
}
1055+
xi_out[index] *= (2*ll3+1.0); //*sign
1056+
}
1057+
} //end omp for
1058+
free(wl_mask);
1059+
free(wigner_sn1);
1060+
if(!same_sn)
1061+
free(wigner_sn2);
1062+
} //end omp parallel
1063+
}
1064+
9971065
void nmt_master_calculator_free(nmt_master_calculator *c)
9981066
{
9991067
int ii, ip, ic;

0 commit comments

Comments
 (0)