Skip to content

Commit

Permalink
Merge pull request #2 from simonsobs/master
Browse files Browse the repository at this point in the history 
updating so-lenspipe
  • Loading branch information
@jaejoonk
jaejoonk authored Sep 26, 2024
2 parents 7514988 + d55cda7 commit ccba364
Show file tree
Hide file tree
Showing 4 changed files with 172 additions and 338 deletions.
12 changes: 12 additions & 0 deletions README.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,18 @@ For a test beyond the imports, you can run
Note that if working on NERSC, you might have to run the scripts on an
interactive node.

Documentation
-------------

This section describes how to provide documentation for the pipeline. You should pip install ``sphinx`` and ``sphinx_rtd_theme`` before proceeding.

* Switch to the `docs` directory and edit the RST files using `reStructured markup <https://sublime-and-sphinx-guide.readthedocs.io/en/latest/index.html>`_.
* Run ``make html`` to build the HTML.
* Open with a browser to view the HTML, e.g. ``firefox _build/html/index.html``
* Submit your edits as a Pull Request.



Contributing
------------

Expand Down
300 changes: 49 additions & 251 deletions solenspipe/bias.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,11 @@
set 2 and 3 should have common phi between corresponding sims, used in MCN1
e.g. implementation:
if set==0 or set==1:
if set==0 or set==1: Corresponding to S and S'
cmb_seed = (icov,set,i)+(0,)
kappa_seed = (icov,set,i)+(1,)
noise_seed = (icov,set,i)+(2,)
elif set==2 or set==3:
elif set==2 or set==3: Corresponds to S_phi and S'_phi of the lensing papers
cmb_seed = (icov,set,i)+(0,)
kappa_seed = (icov,2,i)+(1,)
noise_seed = (icov,set,i)+(2,)
Expand Down Expand Up @@ -626,89 +626,34 @@ def RDN0_analytic(shape,wcs,theory,fwhm,noise_t,noise_p,powdict,estimator,XY,UV,
field_names_alpha=None,field_names_beta=None,skip_filter_field_names=False,
split_estimator=split_estimator)

def mcrdn0_s4(icov, get_kmap, power,phifunc, nsims, qfunc1,get_kmap1=None,get_kmap2=None,get_kmap3=None, qfunc2=None, Xdat=None,Xdat1=None,Xdat2=None,Xdat3=None, use_mpi=True,
def mcrdn0_s4(sim_set,get_kmap, power,phifunc, nsims, qfunc1,get_kmap1=None,get_kmap2=None,get_kmap3=None, qfunc2=None, Xdat=None,Xdat1=None,Xdat2=None,Xdat3=None, use_mpi=True,
verbose=True, skip_rd=False,shear=False,power_mcn0=None):


qa = phifunc
qf1 = qfunc1
qf2=qfunc2


mcn0evals = []
if not(skip_rd):
assert Xdat is not None # Data
if Xdat1 is None:
Xdat1=Xdat
rdn0evals = []

if use_mpi:
comm,rank,my_tasks = mpi.distribute(nsims)
else:
comm,rank,my_tasks = FakeCommunicator(), 0, range(nsims)


for i in my_tasks:
i=i+2
if rank==0 and verbose: print("MCRDN0: Rank %d doing task %d" % (rank,i))
Xs = get_kmap((icov,0,i))
Xs1= get_kmap1((icov,0,i))
Xs2= get_kmap2((icov,0,i))
Xs3= get_kmap3((icov,0,i))


if not(skip_rd):
if shear:
qaXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qbXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf1)
qbXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)
else:
qaXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qbXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf1)
qbXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)

Xsp = get_kmap((icov,0,i+1))
Xsp1 = get_kmap1((icov,0,i+1))
Xsp2 = get_kmap2((icov,0,i+1))
Xsp3 = get_kmap3((icov,0,i+1))

if shear:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs[0],Xsp[1])) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs[0],Xsp[1])) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp[0],Xs[1])) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp[0],Xs[1])) #this is not present
else:
if power_mcn0 is None:
qaXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf1) #split1
qbXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf2) if qf2 is not None else qaXsXsp #split2
qbXspXs = qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf1) #this is not present
mcn0_term = (power(qaXsXsp,qbXsXsp) + power(qaXsXsp,qbXspXs))
else:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs,Xsp)) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs,Xsp)) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp,Xs)) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp,Xs)) #this is not present
mcn0_term = (power_mcn0(qaXsXsp,qbXsXsp) + power_mcn0(qaXsXsp,qbXspXs))

mcn0evals.append(mcn0_term.copy())
if not(skip_rd): rdn0evals.append(rdn0_only_term - mcn0_term)

if not(skip_rd):
avgrdn0 = utils.allgatherv(rdn0evals,comm)
else:
avgrdn0 = None
avgmcn0 = utils.allgatherv(mcn0evals,comm)
return avgrdn0, avgmcn0

"""
This function performs the Realization dependent N0 (RDN0) using combinations of simulations and data for the cross-correlation based estimator. Currently for 4 splits of data.
Parameters:
sim_set: A tuple containing the simulation set and the start index used for the MC. This specifies the simulations used to estimate the RDN0
get_kmap: A function to retrieve a specific map based on the provided parameters.
power: A function to calculate the power spectra
phifunc: function that returns the different phi combinations used to produce a phi_cl without noise bias. Eq.(28)-(30) of 2011.02475
nsims: The number of simulations used to estimate RDN0.
qfunc1: QE wrapper used.
get_kmap1, get_kmap2, get_kmap3: Functions to retrieve the data map splits
qfunc2: An optional QE wrapper if the second phi used to creat clphiphi is different to the first one.
Xdat, Xdat1, Xdat2, Xdat3: Data used in the calculations.
use_mpi: A boolean indicating whether to use MPI for distributing tasks.
verbose: A boolean indicating whether to print verbose messages.
skip_rd: A boolean indicating whether to skip certain operations.
shear: A boolean indicating whether to apply shear in the calculations.
power_mcn0: An optional parameter related to power calculations.
Returns:
A tuple of shape (nsims, 2, Lmax) for the gradient and the curl modes.
def mcrdn0_s41(icov, get_kmap, power,phifunc, nsims, qfunc1,get_kmap1=None,get_kmap2=None,get_kmap3=None, qfunc2=None, Xdat=None,Xdat1=None,Xdat2=None,Xdat3=None, use_mpi=True,
verbose=True, skip_rd=False,shear=False,power_mcn0=None):

"""


i_set,start=sim_set
qa = phifunc
qf1 = qfunc1
qf2=qfunc2
Expand All @@ -728,206 +673,59 @@ def mcrdn0_s41(icov, get_kmap, power,phifunc, nsims, qfunc1,get_kmap1=None,get_k


for i in my_tasks:
i=i+122
i=i+start
if rank==0 and verbose: print("MCRDN0: Rank %d doing task %d" % (rank,i))
Xs = get_kmap((icov,0,i))
Xs1= get_kmap1((icov,0,i))
Xs2= get_kmap2((icov,0,i))
Xs3= get_kmap3((icov,0,i))
Xs = get_kmap((0,i_set,i))
Xs1= get_kmap1((0,i_set,i))
Xs2= get_kmap2((0,i_set,i))
Xs3= get_kmap3((0,i_set,i))


if not(skip_rd):
if shear:
qaXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qaXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf1)
qbXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf1)
qbXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)
else:
qaXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qaXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf1)
qbXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf1)
qbXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)

Xsp = get_kmap((icov,0,i+1))
Xsp1 = get_kmap1((icov,0,i+1))
Xsp2 = get_kmap2((icov,0,i+1))
Xsp3 = get_kmap3((icov,0,i+1))

if shear:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs[0],Xsp[1])) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs[0],Xsp[1])) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp[0],Xs[1])) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp[0],Xs[1])) #this is not present
else:
if power_mcn0 is None:
qaXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf1) #split1
qbXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf2) if qf2 is not None else qaXsXsp #split2
qbXspXs = qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf1) #this is not present
mcn0_term = (power(qaXsXsp,qbXsXsp) + power(qaXsXsp,qbXspXs))
else:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs,Xsp)) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs,Xsp)) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp,Xs)) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp,Xs)) #this is not present
mcn0_term = (power_mcn0(qaXsXsp,qbXsXsp) + power_mcn0(qaXsXsp,qbXspXs))
mcn0evals.append(mcn0_term.copy())
if not(skip_rd): rdn0evals.append(rdn0_only_term - mcn0_term)

if not(skip_rd):
avgrdn0 = utils.allgatherv(rdn0evals,comm)
else:
avgrdn0 = None
avgmcn0 = utils.allgatherv(mcn0evals,comm)
return avgrdn0, avgmcn0


def mcrdn0_s42(icov, get_kmap, power,phifunc, nsims, qfunc1,get_kmap1=None,get_kmap2=None,get_kmap3=None, qfunc2=None, Xdat=None,Xdat1=None,Xdat2=None,Xdat3=None, use_mpi=True,
verbose=True, skip_rd=False,shear=False,power_mcn0=None):


qa = phifunc
qf1 = qfunc1
qf2=qfunc2


mcn0evals = []
if not(skip_rd):
assert Xdat is not None # Data
if Xdat1 is None:
Xdat1=Xdat
rdn0evals = []

if use_mpi:
comm,rank,my_tasks = mpi.distribute(nsims)
else:
comm,rank,my_tasks = FakeCommunicator(), 0, range(nsims)

Xsp = get_kmap((0,i_set,i+1))
Xsp1 = get_kmap1((0,i_set,i+1))
Xsp2 = get_kmap2((0,i_set,i+1))
Xsp3 = get_kmap3((0,i_set,i+1))

for i in my_tasks:
i=i+242
if rank==0 and verbose: print("MCRDN0: Rank %d doing task %d" % (rank,i))
Xs = get_kmap((icov,0,i))
Xs1= get_kmap1((icov,0,i))
Xs2= get_kmap2((icov,0,i))
Xs3= get_kmap3((icov,0,i))

if power_mcn0 is None:

if not(skip_rd):
if shear:
qaXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qbXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf1)
qbXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)
qaXsXsp = plensing.phi_to_kappa(qf1(Xs[0],Xsp[1])) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs[0],Xsp[1])) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp[0],Xs[1])) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp[0],Xs[1])) #this is not present
else:
qaXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qbXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf1)
qbXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)

Xsp = get_kmap((icov,0,i+1))
Xsp1 = get_kmap1((icov,0,i+1))
Xsp2 = get_kmap2((icov,0,i+1))
Xsp3 = get_kmap3((icov,0,i+1))

if shear:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs[0],Xsp[1])) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs[0],Xsp[1])) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp[0],Xs[1])) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp[0],Xs[1])) #this is not present
else:
if power_mcn0 is None:
qaXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf1) #split1
qbXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf2) if qf2 is not None else qaXsXsp #split2
qbXspXs = qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf1) #this is not present
mcn0_term = (power(qaXsXsp,qbXsXsp) + power(qaXsXsp,qbXspXs))
else:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs,Xsp)) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs,Xsp)) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp,Xs)) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp,Xs)) #this is not present
mcn0_term = (power_mcn0(qaXsXsp,qbXsXsp) + power_mcn0(qaXsXsp,qbXspXs))
mcn0evals.append(mcn0_term.copy())
if not(skip_rd): rdn0evals.append(rdn0_only_term - mcn0_term)

if not(skip_rd):
avgrdn0 = utils.allgatherv(rdn0evals,comm)
else:
avgrdn0 = None
avgmcn0 = utils.allgatherv(mcn0evals,comm)
return avgrdn0, avgmcn0


def mcrdn0_s43(icov, get_kmap, power,phifunc, nsims, qfunc1,get_kmap1=None,get_kmap2=None,get_kmap3=None, qfunc2=None, Xdat=None,Xdat1=None,Xdat2=None,Xdat3=None, use_mpi=True,
verbose=True, skip_rd=False,shear=False,power_mcn0=None):


qa = phifunc
qf1 = qfunc1
qf2=qfunc2


mcn0evals = []
if not(skip_rd):
assert Xdat is not None # Data
if Xdat1 is None:
Xdat1=Xdat
rdn0evals = []

if use_mpi:
comm,rank,my_tasks = mpi.distribute(nsims)
else:
comm,rank,my_tasks = FakeCommunicator(), 0, range(nsims)


for i in my_tasks:
i=i+2
if rank==0 and verbose: print("MCRDN0: Rank %d doing task %d" % (rank,i))
Xs = get_kmap((icov,1,i))
Xs1= get_kmap1((icov,1,i))
Xs2= get_kmap2((icov,1,i))
Xs3= get_kmap3((icov,1,i))

mcn0_term = (power(qaXsXsp,qbXsXsp) + power(qaXsXsp,qbXspXs))
else:

if not(skip_rd):
if shear:
qaXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf1)
qbXXs = qa(Xdat[0],Xdat1[0],Xdat2[0],Xdat3[0],Xs[1],Xs1[1],Xs2[1],Xs3[1],qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf1)
qbXsX = qa(Xs[0],Xs1[0],Xs2[0],Xs3[0],Xdat[1],Xdat1[1],Xdat2[1],Xdat3[1],qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)
else:
qaXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf1) #for the two split case, one would need two get_kmaps?, or instead returns an array of maps, [i] for split i
qbXXs = qa(Xdat,Xdat1,Xdat2,Xdat3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qaXXs
qaXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf1)
qbXsX = qa(Xs,Xs1,Xs2,Xs3,Xdat,Xdat1,Xdat2,Xdat3,qf2) if qf2 is not None else qaXsX
rdn0_only_term = power(qaXXs,qbXXs)+ power(qaXXs,qbXsX) + power(qaXsX,qbXXs) \
+ power(qaXsX,qbXsX)

Xsp = get_kmap((icov,1,i+1))
Xsp1 = get_kmap1((icov,1,i+1))
Xsp2 = get_kmap2((icov,1,i+1))
Xsp3 = get_kmap3((icov,1,i+1))

if shear:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs[0],Xsp[1])) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs[0],Xsp[1])) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp[0],Xs[1])) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp[0],Xs[1])) #this is not present
else:
if power_mcn0 is None:
qaXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf1) #split1
qbXsXsp = qa(Xs,Xs1,Xs2,Xs3,Xsp,Xsp1,Xsp2,Xsp3,qf2) if qf2 is not None else qaXsXsp #split2
qbXspXs = qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf2) if qf2 is not None else qa(Xsp,Xsp1,Xsp2,Xsp3,Xs,Xs1,Xs2,Xs3,qf1) #this is not present
mcn0_term = (power(qaXsXsp,qbXsXsp) + power(qaXsXsp,qbXspXs))
qaXsXsp = plensing.phi_to_kappa(qf1(Xs[0],Xsp[1])) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs[0],Xsp[1])) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp[0],Xs[1])) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp[0],Xs[1])) #this is not present
else:
qaXsXsp = plensing.phi_to_kappa(qf1(Xs,Xsp)) #split1
qbXsXsp = plensing.phi_to_kappa(qf2(Xs,Xsp)) if qf2 is not None else qaXsXsp #split2
qbXspXs = plensing.phi_to_kappa(qf2(Xsp,Xs)) if qf2 is not None else plensing.phi_to_kappa(qf1(Xsp,Xs)) #this is not present
mcn0_term = (power_mcn0(qaXsXsp,qbXsXsp) + power_mcn0(qaXsXsp,qbXspXs))
mcn0_term = (power_mcn0(qaXsXsp,qbXsXsp) + power_mcn0(qaXsXsp,qbXspXs))

mcn0evals.append(mcn0_term.copy())
if not(skip_rd): rdn0evals.append(rdn0_only_term - mcn0_term)

Expand Down
Loading

0 comments on commit ccba364

Please sign in to comment.