Rectangle distortion matrix

bin/picca_dmat.py

@@ -47,6 +47,9 @@ def calc_dmat(p):
     parser.add_argument('--nt', type=int, default=50, required=False,
         help='Number of r-transverse bins')
 
+    parser.add_argument('--coef-binning-model', type=int, default=1, required=False,
+        help='Coefficient multiplying np and nt to get finner binning for the model')
+
     parser.add_argument('--z-cut-min', type=float, default=0., required=False,
         help='Use only pairs of forest x object with the mean of the last absorber \
         redshift and the object redshift larger than z-cut-min')
@@ -108,6 +111,8 @@ def calc_dmat(p):
     cf.z_cut_min = args.z_cut_min
     cf.np = args.np
     cf.nt = args.nt
+    cf.npm = args.np*args.coef_binning_model
+    cf.ntm = args.nt*args.coef_binning_model
     cf.nside = args.nside
     cf.zref = args.z_ref
     cf.alpha = args.z_evol
@@ -154,17 +159,29 @@ def calc_dmat(p):
         if not ip in cpu_data:
             cpu_data[ip] = []
         cpu_data[ip].append(p)
-    pool = Pool(processes=args.nproc)
-    dm = pool.map(calc_dmat,sorted(cpu_data.values()))
-    pool.close()
 
+    if args.nproc>1:
+        pool = Pool(processes=args.nproc)
+        dm = pool.map(calc_dmat,sorted(cpu_data.values()))
+        pool.close()
+    elif args.nproc==1:
+        dm = map(calc_dmat,sorted(cpu_data.values()))
+        dm = list(dm)
 
     dm = sp.array(dm)
     wdm =dm[:,0].sum(axis=0)
-    npairs=dm[:,2].sum(axis=0)
-    npairs_used=dm[:,3].sum(axis=0)
+    rp = dm[:,2].sum(axis=0)
+    rt = dm[:,3].sum(axis=0)
+    z = dm[:,4].sum(axis=0)
+    we = dm[:,5].sum(axis=0)
+    npairs = dm[:,6].sum(axis=0)
+    npairs_used = dm[:,7].sum(axis=0)
     dm=dm[:,1].sum(axis=0)
 
+    w = we>0.
+    rp[w] /= we[w]
+    rt[w] /= we[w]
+    z[w] /= we[w]
     w = wdm>0
     dm[w]/=wdm[w,None]
 
@@ -175,11 +192,21 @@ def calc_dmat(p):
         {'name':'RTMAX','value':cf.rt_max,'comment':'Maximum r-transverse [h^-1 Mpc]'},
         {'name':'NP','value':cf.np,'comment':'Number of bins in r-parallel'},
         {'name':'NT','value':cf.nt,'comment':'Number of bins in r-transverse'},
+        {'name':'COEFMOD','value':args.coef_binning_model,'comment':'Coefficient for model binning'},
         {'name':'ZCUTMIN','value':cf.z_cut_min,'comment':'Minimum redshift of pairs'},
         {'name':'ZCUTMAX','value':cf.z_cut_max,'comment':'Maximum redshift of pairs'},
         {'name':'REJ','value':cf.rej,'comment':'Rejection factor'},
         {'name':'NPALL','value':npairs,'comment':'Number of pairs'},
         {'name':'NPUSED','value':npairs_used,'comment':'Number of used pairs'},
     ]
-    out.write([wdm,dm],names=['WDM','DM'],header=head,comment=['Sum of weight','Distortion matrix'],extname='DMAT')
+    out.write([wdm,dm],
+        names=['WDM','DM'],
+        comment=['Sum of weight','Distortion matrix'],
+        units=['',''],
+        header=head,extname='DMAT')
+    out.write([rp,rt,z],
+        names=['RP','RT','Z'],
+        comment=['R-parallel','R-transverse','Redshift'],
+        units=['h^-1 Mpc','h^-1 Mpc','',],
+        extname='ATTRI')
     out.close()

bin/picca_export.py

@@ -91,9 +91,24 @@
     if args.dmat is not None:
         h = fitsio.FITS(args.dmat)
         dm = h[1]['DM'][:]
+        try:
+            dmrp = h[2]['RP'][:]
+            dmrt = h[2]['RT'][:]
+            dmz = h[2]['Z'][:]
+        except IOError:
+            dmrp = rp.copy()
+            dmrt = rt.copy()
+            dmz = z.copy()
+        if dm.shape==(da.size,da.size):
+            dmrp = rp.copy()
+            dmrt = rt.copy()
+            dmz = z.copy()
         h.close()
     else:
         dm = sp.eye(len(da))
+        dmrp = rp.copy()
+        dmrt = rt.copy()
+        dmz = z.copy()
 
     h = fitsio.FITS(args.out,'rw',clobber=True)
     head = [ {'name':'RPMIN','value':rp_min,'comment':'Minimum r-parallel'},
@@ -104,4 +119,6 @@
     ]
     comment = ['R-parallel','R-transverse','Redshift','Correlation','Covariance matrix','Distortion matrix','Number of pairs']
     h.write([rp,rt,z,da,co,dm,nb],names=['RP','RT','Z','DA','CO','DM','NB'],comment=comment,header=head,extname='COR')
+    comment = ['R-parallel model','R-transverse model','Redshift model']
+    h.write([dmrp,dmrt,dmz],names=['DMRP','DMRT','DMZ'],comment=comment,extname='DMATTRI')
     h.close()

bin/picca_export_coadd_zint.py

@@ -42,6 +42,13 @@
 if not args.no_dmat:
     h = fitsio.FITS(args.data[0].replace('cf','dmat'))
     dm = h[1]['DM'][:]*0
+    try:
+        dmrp = sp.zeros(h[2]['RP'][:].size)
+        dmrt = sp.zeros(h[2]['RT'][:].size)
+        dmz = sp.zeros(h[2]['Z'][:].size)
+        nbdm = 0.
+    except IOError:
+        pass
     h.close()
 else:
     dm = sp.eye(nb.shape[0])
@@ -72,6 +79,12 @@
     if not args.no_dmat:
         h = fitsio.FITS(f.replace('cf','dmat'))
         dm += h[1]['DM'][:]*wet_aux[:,None]
+        if 'dmrp' in locals():
+            ## TODO: get the weights
+            dmrp += h[2]['RP'][:]
+            dmrt += h[2]['RT'][:]
+            dmz += h[2]['Z'][:]
+            nbdm += 1.
         h.close()
 
 for p in da:
@@ -91,11 +104,25 @@
 da = (da*we).sum(axis=0)
 da /= wet
 
+
+if 'dmrp' in locals():
+    dmrp /= nbdm
+    dmrt /= nbdm
+    dmz /= nbdm
+elif dmrp.size==rp.size:
+    dmrp = rp.copy()
+    dmrt = rt.copy()
+    dmz = z.copy()
+
 try:
     scipy.linalg.cholesky(co)
 except scipy.linalg.LinAlgError:
     print('WARNING: Matrix is not positive definite')
 
+
 h = fitsio.FITS(args.out,"rw",clobber=True)
-h.write([rp,rt,z,nb,da,dm,co],names=["RP","RT","Z","NB","DA","DM","CO"],header=head,extname='COR')
+comment = ['R-parallel','R-transverse','Redshift','Correlation','Covariance matrix','Distortion matrix','Number of pairs']
+h.write([rp,rt,z,da,co,dm,nb],names=['RP','RT','Z','DA','CO','DM','NB'],comment=comment,header=head,extname='COR')
+comment = ['R-parallel model','R-transverse model','Redshift model']
+h.write([dmrp,dmrt,dmz],names=['DMRP','DMRT','DMZ'],comment=comment,extname='DMATTRI')
 h.close()

bin/picca_metal_dmat.py

@@ -48,6 +48,9 @@ def calc_metal_dmat(abs_igm1,abs_igm2,p):
     parser.add_argument('--nt', type=int, default=50, required=False,
         help='Number of r-transverse bins')
 
+    parser.add_argument('--coef-binning-model', type=int, default=1, required=False,
+        help='Coefficient multiplying np and nt to get finner binning for the model')
+
     parser.add_argument('--z-cut-min', type=float, default=0., required=False,
         help='Use only pairs of forest x object with the mean of the last absorber \
         redshift and the object redshift larger than z-cut-min')
@@ -114,17 +117,16 @@ def calc_metal_dmat(abs_igm1,abs_igm2,p):
     cf.rp_min = args.rp_min
     cf.z_cut_max = args.z_cut_max
     cf.z_cut_min = args.z_cut_min
-    cf.np = args.np
-    cf.nt = args.nt
+    cf.np = args.np*args.coef_binning_model
+    cf.nt = args.nt*args.coef_binning_model
+    cf.npm = args.np*args.coef_binning_model
+    cf.ntm = args.nt*args.coef_binning_model
     cf.nside = args.nside
     cf.zref = args.z_ref
     cf.alpha = args.z_evol
     cf.rej = args.rej
     cf.lambda_abs = constants.absorber_IGM[args.lambda_abs]
     cf.remove_same_half_plate_close_pairs = args.remove_same_half_plate_close_pairs
-    ## use a metal grid equal to the lya grid
-    cf.npm = args.np
-    cf.ntm = args.nt
 
     cf.alpha_abs = {}
     cf.alpha_abs[args.lambda_abs] = cf.alpha
@@ -206,9 +208,15 @@ def calc_metal_dmat(abs_igm1,abs_igm2,p):
             cf.counter.value=0
             f=partial(calc_metal_dmat,abs_igm1,abs_igm2)
             print("")
-            pool = Pool(processes=args.nproc)
-            dm = pool.map(f,sorted(list(cpu_data.values())))
-            pool.close()
+
+            if args.nproc>1:
+                pool = Pool(processes=args.nproc)
+                dm = pool.map(f,sorted(cpu_data.values()))
+                pool.close()
+            elif args.nproc==1:
+                dm = map(f,sorted(cpu_data.values()))
+                dm = list(dm)
+
             dm = sp.array(dm)
             wdm =dm[:,0].sum(axis=0)
             rp = dm[:,2].sum(axis=0)
@@ -241,6 +249,7 @@ def calc_metal_dmat(abs_igm1,abs_igm2,p):
         {'name':'RTMAX','value':cf.rt_max,'comment':'Maximum r-transverse [h^-1 Mpc]'},
         {'name':'NP','value':cf.np,'comment':'Number of bins in r-parallel'},
         {'name':'NT','value':cf.nt,'comment':'Number of bins in r-transverse'},
+        {'name':'COEFMOD','value':args.coef_binning_model,'comment':'Coefficient for model binning'},
         {'name':'ZCUTMIN','value':cf.z_cut_min,'comment':'Minimum redshift of pairs'},
         {'name':'ZCUTMAX','value':cf.z_cut_max,'comment':'Maximum redshift of pairs'},
         {'name':'REJ','value':cf.rej,'comment':'Rejection factor'},

bin/picca_metal_xdmat.py

@@ -44,6 +44,9 @@ def calc_metal_xdmat(abs_igm,p):
     parser.add_argument('--nt', type=int, default=50, required=False,
         help='Number of r-transverse bins')
 
+    parser.add_argument('--coef-binning-model', type=int, default=1, required=False,
+        help='Coefficient multiplying np and nt to get finner binning for the model')
+
     parser.add_argument('--z-min-obj', type=float, default=None, required=False,
         help='Min redshift for object field')
 
@@ -102,17 +105,15 @@ def calc_metal_xdmat(abs_igm,p):
     xcf.rt_max = args.rt_max
     xcf.z_cut_max = args.z_cut_max
     xcf.z_cut_min = args.z_cut_min
-    xcf.np = args.np
-    xcf.nt = args.nt
+    xcf.np = args.np*args.coef_binning_model
+    xcf.nt = args.nt*args.coef_binning_model
+    xcf.npm = args.np*args.coef_binning_model
+    xcf.ntm = args.nt*args.coef_binning_model
     xcf.nside = args.nside
     xcf.zref = args.z_ref
     xcf.lambda_abs = constants.absorber_IGM[args.lambda_abs]
     xcf.rej = args.rej
 
-    ## use a metal grid equal to the lya grid
-    xcf.npm = args.np
-    xcf.ntm = args.nt
-
     cosmo = constants.cosmo(args.fid_Om)
     xcf.cosmo=cosmo
 
@@ -166,9 +167,15 @@ def calc_metal_xdmat(abs_igm,p):
         xcf.counter.value=0
         f=partial(calc_metal_xdmat,abs_igm)
         print("")
-        pool = Pool(processes=args.nproc)
-        dm = pool.map(f,sorted(list(cpu_data.values())))
-        pool.close()
+
+        if args.nproc>1:
+            pool = Pool(processes=args.nproc)
+            dm = pool.map(f,sorted(cpu_data.values()))
+            pool.close()
+        elif args.nproc==1:
+            dm = map(f,sorted(cpu_data.values()))
+            dm = list(dm)
+
         dm = sp.array(dm)
         wdm =dm[:,0].sum(axis=0)
         rp = dm[:,2].sum(axis=0)
@@ -201,6 +208,7 @@ def calc_metal_xdmat(abs_igm,p):
         {'name':'RTMAX','value':xcf.rt_max,'comment':'Maximum r-transverse [h^-1 Mpc]'},
         {'name':'NP','value':xcf.np,'comment':'Number of bins in r-parallel'},
         {'name':'NT','value':xcf.nt,'comment':'Number of bins in r-transverse'},
+        {'name':'COEFMOD','value':args.coef_binning_model,'comment':'Coefficient for model binning'},
         {'name':'ZCUTMIN','value':xcf.z_cut_min,'comment':'Minimum redshift of pairs'},
         {'name':'ZCUTMAX','value':xcf.z_cut_max,'comment':'Maximum redshift of pairs'},
         {'name':'REJ','value':xcf.rej,'comment':'Rejection factor'},

bin/picca_xdmat.py

@@ -43,6 +43,9 @@ def calc_dmat(p):
     parser.add_argument('--nt', type=int, default=50, required=False,
         help='Number of r-transverse bins')
 
+    parser.add_argument('--coef-binning-model', type=int, default=1, required=False,
+        help='Coefficient multiplying np and nt to get finner binning for the model')
+
     parser.add_argument('--z-min-obj', type=float, default=None, required=False,
         help='Min redshift for object field')
 
@@ -99,6 +102,8 @@ def calc_dmat(p):
     xcf.z_cut_min = args.z_cut_min
     xcf.np = args.np
     xcf.nt = args.nt
+    xcf.npm = args.np*args.coef_binning_model
+    xcf.ntm = args.nt*args.coef_binning_model
     xcf.nside = args.nside
     xcf.zref = args.z_ref
     xcf.alpha = args.z_evol_del
@@ -151,30 +156,52 @@ def calc_dmat(p):
             cpu_data[ip] = []
         cpu_data[ip].append(p)
 
-    pool = Pool(processes=args.nproc)
-    dm = pool.map(calc_dmat,sorted(list(cpu_data.values())))
-    pool.close()
+    if args.nproc>1:
+        pool = Pool(processes=args.nproc)
+        dm = pool.map(calc_dmat,sorted(list(cpu_data.values())))
+        pool.close()
+    elif args.nproc==1:
+        dm = map(calc_dmat,sorted(list(cpu_data.values())))
+        dm = list(dm)
+
     dm = sp.array(dm)
     wdm =dm[:,0].sum(axis=0)
-    npairs=dm[:,2].sum(axis=0)
-    npairs_used=dm[:,3].sum(axis=0)
+    rp = dm[:,2].sum(axis=0)
+    rt = dm[:,3].sum(axis=0)
+    z = dm[:,4].sum(axis=0)
+    we = dm[:,5].sum(axis=0)
+    npairs = dm[:,6].sum(axis=0)
+    npairs_used = dm[:,7].sum(axis=0)
     dm=dm[:,1].sum(axis=0)
 
+    w = we>0.
+    rp[w] /= we[w]
+    rt[w] /= we[w]
+    z[w] /= we[w]
     w = wdm>0.
     dm[w,:] /= wdm[w,None]
 
-
     out = fitsio.FITS(args.out,'rw',clobber=True)
     head = [ {'name':'RPMIN','value':xcf.rp_min,'comment':'Minimum r-parallel [h^-1 Mpc]'},
         {'name':'RPMAX','value':xcf.rp_max,'comment':'Maximum r-parallel [h^-1 Mpc]'},
         {'name':'RTMAX','value':xcf.rt_max,'comment':'Maximum r-transverse [h^-1 Mpc]'},
         {'name':'NP','value':xcf.np,'comment':'Number of bins in r-parallel'},
         {'name':'NT','value':xcf.nt,'comment':'Number of bins in r-transverse'},
+        {'name':'COEFMOD','value':args.coef_binning_model,'comment':'Coefficient for model binning'},
         {'name':'ZCUTMIN','value':xcf.z_cut_min,'comment':'Minimum redshift of pairs'},
         {'name':'ZCUTMAX','value':xcf.z_cut_max,'comment':'Maximum redshift of pairs'},
         {'name':'REJ','value':xcf.rej,'comment':'Rejection factor'},
         {'name':'NPALL','value':npairs,'comment':'Number of pairs'},
         {'name':'NPUSED','value':npairs_used,'comment':'Number of used pairs'},
     ]
-    out.write([wdm,dm],names=['WDM','DM'],header=head,comment=['Sum of weight','Distortion matrix'],extname='DMAT')
+    out.write([wdm,dm],
+        names=['WDM','DM'],
+        comment=['Sum of weight','Distortion matrix'],
+        units=['',''],
+        header=head,extname='DMAT')
+    out.write([rp,rt,z],
+        names=['RP','RT','Z'],
+        comment=['R-parallel','R-transverse','Redshift'],
+        units=['h^-1 Mpc','h^-1 Mpc','',],
+        extname='ATTRI')
     out.close()