<jupyter_start><jupyter_text>Fitting cosmological parameters with the Core Cosmological Library (CCL)This notebook shows how to fit cosmological parameters using CCL and measurements of cosmic shear correlation functions ($\xi_\pm$). Currently it uses analytical forms for the photometric redshift count histograms, $n(z)$ and simulates data based on the correlation functions computed by CCL. Replacing [$n(z)$](nz) and [$\xi_\pm^{ij}$](datasets) with real data should be straight forward.<jupyter_code>%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines){
    return false;
}
import pylab as plt
from matplotlib import gridspec
import numpy as np
import ipywidgets as widgets
from ipywidgets import interact
import itertools
import pyccl as ccl
import os
import scipy.stats as st
from iminuit import Minuit
from iminuit.util import describe, make_func_code
from iminuit.cost import Cost

%matplotlib inline<jupyter_output><empty_output><jupyter_text>`Cosmo` class by <NAME>. Reparametrizes cosmological parameters to be non-degenerate for cosmic shear measurements.<jupyter_code>class Cosmo:
    def __init__(self,Omega_ch2=0.1, Omega_bh2=0.023, AS=None, S8=None,
                 Omega_nu_h2=0., H0=70, ns=0.97, w0=-1, alpha=0.45, 
                 matter_power_spectrum = 'halofit'):
        
        self._matter_power_spectrum =  matter_power_spectrum
        self.update_cosmology(Omega_ch2=Omega_ch2, Omega_bh2=Omega_bh2,
                              AS=AS, S8=S8, alpha=alpha,
                              Omega_nu_h2=Omega_nu_h2,
                              H0=H0, ns=ns, w0=w0)
        
    def update_cosmology(self, Omega_ch2=0.1, Omega_bh2=0.023, AS=None, S8=None,
                         Omega_nu_h2=1e-3, H0=70, ns=0.97, w0=-1, alpha=0.45):

        self.Omega_ch2 = Omega_ch2
        self.Omega_bh2 = Omega_bh2
        self.H0 = H0
        self._h = self.H0 / 100.

        self._Omega_b = self.Omega_bh2 / self._h**2
        self._Omega_c = self.Omega_ch2 / self._h**2
        self._Omega_m = self._Omega_b + self._Omega_c
        self.Omega_nu_h2 = Omega_nu_h2
        self._m_nu = (self.Omega_nu_h2 / self._h**2) * 93.14

        if S8 is None and AS is None:
            raise ValueError('S8 or AS should be given')
        if S8 is not None and AS is not None:
            raise ValueError('Just S8 or AS should be given')

        if S8 is not None:
            self.AS = None
            self._A_s = None
            self._sigma8 = S8 * (1./ (self._Omega_m/0.3)**alpha)
        if AS is not None:
            self.AS = AS
            self._A_s = np.exp(self.AS) * 1e-10
            self._sigma8 = None

        self.n_s = ns
        self.w0 = w0

        self.cosmology = ccl.Cosmology(Omega_c=self._Omega_c, Omega_b=self._Omega_b,
                                       h=self._h, n_s=self.n_s, sigma8=self._sigma8, A_s=self._A_s,
                                       w0=self.w0, m_nu=self._m_nu,
                                       matter_power_spectrum=self._matter_power_spectrum)

C = Cosmo(Omega_ch2=0.1, Omega_bh2=0.023, AS=4.,
          Omega_nu_h2=1e-3, H0=70, ns=0.97, w0=-1,
          matter_power_spectrum = 'halofit')
C.cosmology
@interact(Om = widgets.FloatSlider(value=0.9, min=0.01, max=0.99, step=0.01, description='$\Omega_m$:',
          disabled=False,
          continuous_update=False,
          orientation='horizontal'),
          w = widgets.FloatSlider(value=-1., min=-2, max=-0.3, step=0.01, description='$w$:',
          disabled=False,
          continuous_update=False,
          orientation='horizontal',
          readout=True,
          readout_format='.2f'))
def distance_modulus(Om, w):
          
    C = Cosmo(Omega_ch2=Om*0.7*0.7, Omega_bh2=0.0, AS=4.,
              Omega_nu_h2=0., H0=70, ns=0.97, w0=w,
              matter_power_spectrum = 'halofit')
    z = np.linspace(0.01, 1.2, 200)
    plt.figure(figsize=(10,6))
    plt.plot(z, ccl.distance_modulus(C.cosmology, 1/(1+z)))
    plt.xlabel('z',fontsize=18)
    plt.ylabel('$\mu$', fontsize=18)
    plt.ylim(32,45)
    plt.xlim(0.01,1.2)
    #plt.xscale('log')<jupyter_output><empty_output><jupyter_text>All weak lensing tracers use the following analytic Gaussian redshift bins. Note `redshifts` and `nz` are global variables.<jupyter_code>key = '_hamana'

redshifts = []
nz = []

# for i in range(4):
#     file_bin = np.loadtxt(os.path.join('', '../Analysis/shear_subaru/shear_subaru/data/photo-z/bin%i'%(i+1))+key+'.dat', comments='#')
#     redshifts.append(file_bin[:,0])
#     nz.append(file_bin[:,1])
    
for i in range(4):
    z_range = np.linspace(0, 2.6, 2600)
    nz_i = st.norm.pdf(z_range, (i+1)*0.4, 0.2)
    redshifts.append(z_range)
    nz.append(nz_i)

redshifts = np.array(redshifts)
nz_distros = np.array(nz)

@interact(dz1 = widgets.FloatSlider(value=0., min=-0.1, max=0.1, step=0.01, description='$\Delta z_1$:',
          disabled=False,
          continuous_update=False,
          orientation='horizontal'),
          dz2 = widgets.FloatSlider(value=0., min=-0.1, max=0.1, step=0.01, description='$\Delta z_2$:',
          disabled=False,
          continuous_update=False,
          orientation='horizontal'),
          dz3 = widgets.FloatSlider(value=0., min=-0.1, max=0.1, step=0.01, description='$\Delta z_3$:',
          disabled=False,
          continuous_update=False,
          orientation='horizontal'),
          dz4 = widgets.FloatSlider(value=0., min=-0.1, max=0.1, step=0.01, description='$\Delta z_4$:',
          disabled=False,
          continuous_update=False,
          orientation='horizontal'))

def load_photo_z(dz1, dz2, dz3, dz4):
    
    DZ = [dz1, dz2, dz3, dz4]

 
    C = ['k', 'b', 'y', 'r']
    plt.figure(figsize=(12,4))
    for i in range(4):
        plt.plot(redshifts[i]-DZ[i], nz[i], C[i], lw=3)

    plt.plot([0,2.6], [0,0], 'k--')
    plt.xlabel('z', fontsize=18)
    plt.ylabel('N(z)', fontsize=18)
    plt.xlim(0,2.6)
    plt.ylim(-0.1,4.2)<jupyter_output><empty_output><jupyter_text>Set up weak lensing tracers for each redshift bin as a demonstration. These instances don't actually get used for fitting.<jupyter_code>def intrinsic_al(redshift, A0=1, eta=1, z0=0.62):
    AI = A0 * ((1+redshift) / (1+z0))**eta
    return AI

WL_bin1 = ccl.WeakLensingTracer(C.cosmology,
                                dndz=(redshifts[0], nz[0]),
                                has_shear=True,
                                ia_bias=(redshifts[0], intrinsic_al(redshifts[0])))

WL_bin2 = ccl.WeakLensingTracer(C.cosmology,
                                dndz=(redshifts[1], nz[1]),
                                has_shear=True,
                                ia_bias=(redshifts[1], intrinsic_al(redshifts[1])))

WL_bin3 = ccl.WeakLensingTracer(C.cosmology,
                                dndz=(redshifts[2], nz[2]),
                                has_shear=True,
                                ia_bias=(redshifts[2], intrinsic_al(redshifts[2])))

WL_bin4 = ccl.WeakLensingTracer(C.cosmology,
                                dndz=(redshifts[3], nz[3]),
                                has_shear=True,
                                ia_bias=(redshifts[3], intrinsic_al(redshifts[3])))<jupyter_output>/home/yve/anaconda3/lib/python3.7/site-packages/pyccl/core.py:678: CCLWarning: CCL does not properly compute the linear growth rate in cosmological models with massive neutrinos!
  category=CCLWarning)<jupyter_text>**Main model function**<jupyter_code>def xi_pm_ij(theta, Om_ch2, AS, bin_i=1, bin_j=1, sign='+'):
    """
    Currently can fit for Om_ch2, AS
    bin_i, bin_j, and sign must be passed through the cost function: is
    there a better way to generalize this?
    """
    
    #set up cosmology to be fit
    #currently only fits for Omega_ch2 and AS
    C = Cosmo(Omega_ch2=Om_ch2, Omega_bh2=0.023, AS=AS,
          Omega_nu_h2=1e-3, H0=70, ns=0.97, w0=-1,
          matter_power_spectrum = 'halofit')
        
    #set up redshift and n(z) bins from (global variable) data
    bin_i_idx = int(bin_i - 1)
    bin_j_idx = int(bin_j - 1)
    
    z_i = redshifts[bin_i_idx]
    nz_i = nz_distros[bin_i_idx]
    
    z_j = redshifts[bin_j_idx]
    nz_j = nz_distros[bin_j_idx]
    
    AI_i = intrinsic_al(z_i, A0=1, eta=1, z0=0.62)
    AI_j = intrinsic_al(z_j, A0=1, eta=1, z0=0.62)

    #calculate WL bins
    WL_bin_i = ccl.WeakLensingTracer(C.cosmology,
                                    dndz=(z_i, nz_i),
                                    has_shear=True,
                                    ia_bias=(z_i, AI_i))
    
    WL_bin_j = ccl.WeakLensingTracer(C.cosmology,
                                    dndz=(z_j, nz_j),
                                    has_shear=True,
                                    ia_bias=(z_j, AI_j))
    
    #Calculate angular power spectrum
    ell = np.arange(1, 2000)
    Cl_ij = ccl.angular_cl(C.cosmology, WL_bin_i, WL_bin_j, ell)

    #calculate correlation function xi_plus_ij or xi_minus_ij
    if sign is '+':
        type_str = 'GG+'
    else:
        type_str = 'GG-'
    
    xi_pm_ij = ccl.correlation(C.cosmology, ell, Cl_ij, theta, type=type_str)
        
    return xi_pm_ij<jupyter_output><empty_output><jupyter_text>**Preliminary cost function for data with no errors.**\We will see this doesn't behave very well...<jupyter_code>#adapted class from iminuit tutorial
class LeastSquaresNoErr(Cost):
    """
    Generic least-squares cost function without error.
    """

    errordef = Minuit.LEAST_SQUARES # for Minuit to compute errors correctly
    
    def __init__(self, model, x, y, bin_i=None, bin_j=None, sign=None):
        super().__init__(describe(model)[1:3], len(x), verbose=1)
        self.model = model  # model predicts y for given x
        self.x = np.asarray(x)
        self.y = np.asarray(y)
        self.bin_i = bin_i
        self.bin_j = bin_j
        self.sign = sign        
        
    def _call(self, *par):  # we accept a variable number of model parameters
        ym = self.model(self.x, bin_i = self.bin_i, bin_j = self.bin_j, sign=self.sign, *par[0])
        return np.sum((self.y - ym) ** 2)<jupyter_output><empty_output><jupyter_text>Produce data for Xi_44_plus for testing model and cost functions<jupyter_code>#get samples from xip_44
ell = np.arange(1, 2000)
theta_sample = np.logspace(-2, np.log10(2), 20)

Cl_44 = ccl.angular_cl(C.cosmology, WL_bin4, WL_bin4, ell)
xip_44_sample = ccl.correlation(C.cosmology, ell, Cl_44, theta_sample, type='GG+')
print(theta_sample)
print(xip_44_sample)

#add noise
noise = st.norm.rvs(0, 5e-5, size=len(xip_44_sample), random_state=1369)
noisy_44_sample = xip_44_sample + noise

print(noisy_44_sample)


plt.figure()
plt.scatter(60*theta_sample, xip_44_sample, label='44')
plt.scatter(60*theta_sample, noisy_44_sample, label='44+noise')
plt.xscale('log')
plt.xlabel('$\\theta$', fontsize=18)
plt.ylabel('$\\xi_{+}^{44}$', fontsize=18)
plt.ylim(-0.0001, np.max(xip_44_sample)+0.1*np.max(xip_44_sample))
plt.legend()
#Testing Least-Squares with no error on Xi_plus_44 correlation "data" with no errors/noise.
lsq_no_err = LeastSquaresNoErr(xi_pm_ij, theta_sample, xip_44_sample, 4, 4, '+')

#input cosmology for reference
#C = Cosmo(Omega_ch2=0.1, Omega_bh2=0.023, AS=4.,
#          Omega_nu_h2=0, H0=70, ns=0.97, w0=-1,
#          matter_power_spectrum = 'halofit')

m = Minuit(lsq_no_err, Om_ch2=0.11, AS=4.1)
m.params

m.limits['Om_ch2'] = (0.03, 0.7)
m.limits['AS'] = (1.5, 6)

m.migrad()
print(m.valid, m.accurate)
m.draw_profile('Om_ch2', bound=(0.03, 0.6), size=20)
plt.yscale('log')
m.draw_profile('AS', bound=(1.5, 6), size=20)
plt.yscale('log')<jupyter_output>/home/yve/anaconda3/lib/python3.7/site-packages/pyccl/core.py:678: CCLWarning: CCL does not properly compute the linear growth rate in cosmological models with massive neutrinos!
  category=CCLWarning)<jupyter_text>Plotting contours for this set fails because `draw_contour` expects there will be contours for cost function values of 1,2,3, and 4, and as we see above, the cost function has values << 1.<jupyter_code>#plot Om_ch2, AS contour plot
# x, y = 'Om_ch2', 'AS'

# m.draw_contour(x, y, size=10, bound=[[0.03,0.2],[3.,5.]])
# plt.yscale('log')<jupyter_output><empty_output><jupyter_text>Create datasets We now create dictionaries for datasets for all auto- and cross-correlation functions.\To summarize all dictionaries created in next 3 cells:\`data_dict` = Exact data calculated from Cosmology defined at top of notebook.\`sigma_arr_dict` = Exact 10% of each data point in `data_dict`. Used as sigma to produce noise for following dictionaries.\`seed_noisy_data_dict` = data where all datapoints for any function (e.g. Xi_14_+) have ~10% Gaussian errors **using the same seed for all functions.**\`random_noisy_data_dict` = data where where all datapoints for any function have ~10% **random** Gaussian errors.<jupyter_code>#create dictionary of N data points from each xi_plus and xi_minus function
N = 20
ell = np.arange(1, 2000)
theta_sample = np.logspace(-2, np.log10(2), N)

data_dict = {}

for i in range(4):
    for j in range(i,4,1):
        pair_str = str(i+1)+str(j+1)
        xip_str = 'xip_'+pair_str
        xim_str = 'xim_'+pair_str
        
        xi_p_ij = xi_pm_ij(theta_sample, 0.1, 4, i+1, j+1, '+')
        xi_m_ij = xi_pm_ij(theta_sample, 0.1, 4, i+1, j+1, '-')
        
        data_dict[xip_str] = xi_p_ij
        data_dict[xim_str] = xi_m_ij

data_dict
#create noisy data with ~10% errors using fixed random seed for all correlation functions
import copy

seed_noisy_data_dict = {}
sigma_arr_dict = copy.deepcopy(data_dict) #dictionary of iterables: must use copy.deepcopy()

for key in data_dict:
    sigma_arr_dict[key] /= 10 #ten percent errors for each data point

    noise = st.norm.rvs(0, sigma_arr_dict[key], random_state=1369)
    xi_ij = data_dict[key]
    
    seed_noisy_data_dict[key] = xi_ij + noise
    
seed_noisy_data_dict
#make correlation function dataset with random Gaussian ~10% errors
random_noisy_data_dict = {}
for key in data_dict:
    random_noisy_data_dict[key] = st.norm.rvs(data_dict[key], sigma_arr_dict[key])
    
random_noisy_data_dict<jupyter_output><empty_output><jupyter_text>Joint fit all 1-4 bin combinations (20 correlation fctns) Uses data with no errors. Expensive to run!\I lost the output for this one; I leave running it as an exercise to the reader... :P<jupyter_code>lsqne_tot = None
for i in range(4):
    for j in range(i,4,1):
        pair_str = str(i+1)+str(j+1)
        xip_str = 'xip_'+pair_str
        xim_str = 'xim_'+pair_str
        lsqne_ij = (LeastSquaresNoErr(xi_pm_ij, theta_sample, data_dict[xip_str], i+1, j+1, '+')
                   +LeastSquaresNoErr(xi_pm_ij, theta_sample, data_dict[xim_str], i+1, j+1, '-'))
        
        if lsqne_tot is not None:
            lsqne_tot += lsqne_ij
        else:
            lsqne_tot = lsqne_ij

m2 = Minuit(lsqne_tot, Om_ch2=0.11, AS=4.1)

m2.limits['Om_ch2'] = (0.03, 0.7)
m2.limits['AS'] = (1.5, 6)

m2.migrad()<jupyter_output><empty_output><jupyter_text>Custom Least Squares cost function: Non-linear least squares tailed to CCL shear-shear correlation function (`xi_ij_pm`).\*NB:* this currently only fits exactly 2 parameters. For $k$ parameters, use `describe(model)[1:k+1]` in `__init__`.<jupyter_code>class LeastSquares(Cost):
    """
    Generic least-squares cost function tailored for CCL shear-shear correlation function.
    """

    errordef = Minuit.LEAST_SQUARES # for Minuit to compute errors correctly
    
    def __init__(self, model, x, y, yerr):
        super().__init__(describe(model)[1:], len(x), verbose=1)
        self.model = model  # model predicts y for given x
        self.x = np.asarray(x)
        self.y = np.asarray(y)
        if len(y)==1:
            self.yerr = yerr * np.ones_like(y)
        else:
            self.yerr = yerr
    
    #overload __init__ for CCL model function
    def __init__(self, model, x, y, yerr, bin_i=None, bin_j=None, sign=None):
        super().__init__(describe(model)[1:3], len(x), verbose=1) #this will need to be changed 
                                                                  #if adding more fitting parameters
        self.model = model
        self.x = np.asarray(x)
        self.y = np.asarray(y)
        if len(y)==1:
            self.yerr = yerr * np.ones_like(y)
        else:
            self.yerr = np.asarray(yerr)       
        self.bin_i = bin_i
        self.bin_j = bin_j
        self.sign = sign
        
    def _call(self, *par):  # we accept a variable number of model parameters
        ym = self.model(self.x, bin_i = self.bin_i, bin_j = self.bin_j, sign=self.sign, *par[0])
        #ym = self.model(self.x, *par[0]) #uncomment to test linear data below
        return np.sum(((self.y - ym) / self.yerr) ** 2)
#test with linear data
def linear(x, a, b):
    return a*x + b

x = np.linspace(0,10,10)
a = 2
b = 3
y = st.norm.rvs(a*x+b, 2, random_state=720)

plt.plot(x, linear(x, a, b))
plt.errorbar(x, y, 2, fmt='o')
lsq_lin = LeastSquares(linear, x, y, 2)
mlin = Minuit(lsq_lin, a=1, b=1)

mlin.migrad()
plt.figure(figsize=(10,6))
plt.plot(x, linear(x, mlin.values['a'], mlin.values['b']), label='migrad fit')
plt.plot(x, linear(x, a, b), label='Mean for generating model')
plt.errorbar(x, y, 2, fmt='o', label='generated data')
plt.legend()
plt.figure(figsize=(10,6))
mlin.draw_profile('a', subtract_min=True)
plt.axvline(1.63, color='r', label=r'$\Delta\chi^2 = 1$: 1.63-1.847')
plt.axvline(1.847, color='r')

plt.plot(np.linspace(1.34,2.13,2), np.ones(2), linestyle='--')
plt.legend()
plt.figure(figsize=(10,6))
mlin.draw_profile('b', subtract_min=True)
plt.axvline(4.41, color='r', label=r'$\Delta\chi^2 = 1$: 4.41-5.70')
plt.axvline(5.70, color='r')
plt.plot(np.linspace(2.7,7.41,2), np.ones(2), linestyle='--')
plt.legend()
mlin.draw_contour('a', 'b', size=10)<jupyter_output>(1.342113395100168, 2.7091731582041256) -> 60.27474905066536
(1.342113395100168, 3.2316209937997424) -> 49.63894899411339
(1.342113395100168, 3.7540688293953597) -> 40.36790764215412
(1.342113395100168, 4.276516664990977) -> 32.46162499478759
(1.342113395100168, 4.798964500586594) -> 25.92010105201377
(1.342113395100168, 5.321412336182211) -> 20.743335813832694
(1.342113395100168, 5.843860171777828) -> 16.931329280244334
(1.342113395100168, 6.366308007373444) -> 14.484081451248699
(1.342113395100168, 6.888755842969061) -> 13.401592326845787
(1.342113395100168, 7.411203678564679) -> 13.683861907035594
(1.4301904965271188, 2.7091731582041256) -> 49.63894899323629
(1.4301904965271188, 3.2316209937997424) -> 40.153541211835474
(1.4301904965271188, 3.7540688293953597) -> 32.03289213502735
(1.4301904965271188, 4.276516664990977) -> 25.27700176281197
(1.4301904965271188, 4.798964500586594) -> 19.885870095189308
(1.4301904965271188, 5.321412336182211) -> 15.859497132159381
(1.4301904965271188,[...]<jupyter_text>Comparing with the LeastSquares function packaged with iminuit.\Gets same results as above:<jupyter_code>from iminuit.cost import LeastSquares as minuitLS

minuit_lsq_lin = minuitLS(x, y, 2, linear, verbose=1)
minuit_mlin = Minuit(minuit_lsq_lin, a=1, b=1)

minuit_mlin.migrad()
plt.figure(figsize=(10,6))
minuit_mlin.draw_profile('a', subtract_min=True)
plt.axvline(1.63, color='r', label='1.63')
plt.axvline(1.847, color='r', label='1.847')
plt.plot(np.linspace(1.34,2.13,2), np.ones(2), linestyle='--')
plt.title('Using minuit\'s Least-Squares - same result')
plt.legend()
plt.figure(figsize=(10,6))
minuit_mlin.draw_profile('b', subtract_min=True)
plt.axvline(4.41, color='r', label='4.41')
plt.axvline(5.70, color='r', label='5.70')
plt.plot(np.linspace(2.7,7.41,2), np.ones(2), linestyle='--')
plt.legend()<jupyter_output>(1.7384603515213157, 2.7091731739716027) -> 23.161123407288713
(1.7384603515213157, 2.756668431434482) -> 22.608452534329455
(1.7384603515213157, 2.804163688897361) -> 22.067060658777496
(1.7384603515213157, 2.8516589463602395) -> 21.536947780632897
(1.7384603515213157, 2.8991542038231186) -> 21.01811389989559
(1.7384603515213157, 2.9466494612859977) -> 20.51055901656563
(1.7384603515213157, 2.9941447187488768) -> 20.014283130642976
(1.7384603515213157, 3.041639976211756) -> 19.529286242127675
(1.7384603515213157, 3.0891352336746345) -> 19.05556835101968
(1.7384603515213157, 3.1366304911375136) -> 18.593129457319016
(1.7384603515213157, 3.1841257486003927) -> 18.141969561025675
(1.7384603515213157, 3.2316210060632717) -> 17.702088662139666
(1.7384603515213157, 3.2791162635261504) -> 17.273486760660983
(1.7384603515213157, 3.3266115209890295) -> 16.856163856589614
(1.7384603515213157, 3.3741067784519085) -> 16.45011994992558
(1.7384603515213157, 3.4216020359147876) -> 16.05535504066887
[...]<jupyter_text>Set up test for non-linear least squares with CCL model<jupyter_code>test_str = 'xip_44'
plt.figure(figsize=(10,6))
np.logspace(-2, np.log10(2), 20)
theta = np.logspace(-3, np.log10(2), 100)
xip_44 = xi_pm_ij(theta, 0.1, 4., 4, 4, '+')
plt.plot(theta, xip_44, label='$\\xi_{+}^{44}$')
plt.errorbar(theta_sample, seed_noisy_data_dict[test_str], sigma_arr_dict[test_str], 
             fmt='o', label='sample w/ 10% seeded errors')
plt.errorbar(theta_sample, random_noisy_data_dict[test_str], sigma_arr_dict[test_str], 
             fmt='o', label='sample w/ 10% random errors')

plt.ylim((5e-6,0.002))
plt.xscale('log')
plt.yscale('log')
plt.legend()<jupyter_output>/home/yve/anaconda3/lib/python3.7/site-packages/pyccl/core.py:678: CCLWarning: CCL does not properly compute the linear growth rate in cosmological models with massive neutrinos!
  category=CCLWarning)<jupyter_text>Fitting with fixed seed data First test with seeded data (orange data points).<jupyter_code>#Testing custom LeastSquares with Xi_44_plus auto-correlation data with seeded Gaussian noise
lsq_p_44 = LeastSquares(xi_pm_ij, theta_sample, seed_noisy_data_dict[test_str], 
                        sigma_arr_dict[test_str], 4, 4, '+')

m3 = Minuit(lsq_p_44, Om_ch2=0.11, AS=4.1)

m3.limits['Om_ch2'] = (0.03, 0.7)
m3.limits['AS'] = (1.5, 6)

m3.migrad()
m3.draw_profile('Om_ch2', bound=5)
# plt.yscale('log')
m3.draw_profile('AS', bound=5, size=25)
# plt.yscale('log')
m3.draw_contour('Om_ch2', 'AS', size=10, bound=[[0.08,0.1],[4.05,4.25]])<jupyter_output>/home/yve/anaconda3/lib/python3.7/site-packages/pyccl/core.py:678: CCLWarning: CCL does not properly compute the linear growth rate in cosmological models with massive neutrinos!
  category=CCLWarning)<jupyter_text>We can also use the `minos` optimization algorithm to estimate parameters. This is more computationally expensive than `migrad` (see iminuit docs).<jupyter_code>m3.minos()
om_results = m3.mnprofile('Om_ch2', size=10, bound=(0.08, 0.1))
omch2_arr = om_results[0]
fval_arr = om_results[1]
delta_fval_arr = fval_arr - np.min(fval_arr)

plt.plot(omch2_arr, delta_fval_arr)
m3.mnprofile('AS', size=10, bound=(4.0, 4.2))
as_arr = np.array([4., 4.02222222, 4.04444444, 4.06666667, 4.08888889,
        4.11111111, 4.13333333, 4.15555556, 4.17777778, 4.2])
fval_2_arr = np.array([7.85418292, 7.32111812, 6.86757683, 6.4945691 , 6.20338035,
        5.99511007, 5.8724382 , 5.83099018, 5.87851035, 6.01368415])
fval_2_min_arr = fval_2_arr - np.min(fval_2_arr)

plt.plot(as_arr, fval_2_min_arr)
m3.mncontour('Om_ch2', 'AS', size=10)
m3_contour = np.array([[0.08283136, 4.31018723],
       [0.08543363, 4.23016267],
       [0.08891746, 4.15017163],
       [0.09293852, 4.07020434],
       [0.098267  , 3.99031351],
       [0.09841201, 3.9930902 ],
       [0.0953999 , 4.07315942],
       [0.09166628, 4.15319825],
       [0.08771533, 4.23322685],
       [0.08295866, 4.3132116 ]])

plt.plot(contour[:,0], contour[:,1])<jupyter_output><empty_output><jupyter_text>Fitting with perfect data w/ ~10% errors Testing migrad and minos with perfect data (`data_dict`) w/ ~10% errors.<jupyter_code>lsq_p_44_v3 = LeastSquares(xi_pm_ij, theta_sample, data_dict['xip_44'], 
                           sigma_arr_dict['xip_44'], 4, 4, '+')

m4 = Minuit(lsq_p_44_v3, Om_ch2=0.11, AS=4.1)

m4.limits['Om_ch2'] = (0.03, 0.7)
m4.limits['AS'] = (1.5, 6)

m4.migrad()
m4.draw_profile('Om_ch2', bound=2)
m4.draw_profile('AS', bound=2)
m4.draw_contour('Om_ch2', 'AS', size=10)
m4.minos()<jupyter_output>(0.10586403372825617, 3.8917645468828477) -> 0.9812532955557385
(0.10586403372825617, 3.8928605378957566) -> 0.9656424183496629
(0.10586403372825617, 3.890668522097884) -> 1.0032572800492536
(0.10586403372825617, 3.8918741475081893) -> 0.9794038504639717
(0.10586403372825617, 3.8916549459197856) -> 0.9831667336558652
(0.10586403372825617, 3.8949581309400223) -> 0.9536953201221614
(0.10586403372825617, 3.8949713536777093) -> 0.9536949844666753
(0.10586403372825617, 3.8949449081973073) -> 0.953696597397748
(0.10600459148157547, 3.8949581309400223) -> 1.0141676268958202
(0.10600459148157547, 3.8951513534472717) -> 1.0168150999009193
(0.10600459148157547, 3.894764907359273) -> 1.011722779143164
(0.10600459148157547, 3.894977453239078) -> 1.0144232496926529
(0.10600459148157547, 3.894938808630231) -> 1.0139139362385072
(0.10600459148157547, 3.892497767769216) -> 0.9981359060022179
(0.10600459148157547, 3.8925113373869378) -> 0.9981346793822359
(0.10600459148157547, 3.892484198146291) -> 0.9[...]<jupyter_text>Fitting to $\xi_+^{44}$  and  $\xi_-^{44}$<jupyter_code>#test adding two cost functions
xip_str = 'xip_44'
xim_str = 'xim_44'
lsq_pm44 = (LeastSquares(xi_pm_ij, theta_sample, random_noisy_data_dict[xip_str], 
                       sigma_arr_dict[xip_str], 4, 4, '+')
          +LeastSquares(xi_pm_ij, theta_sample, random_noisy_data_dict[xim_str],
                        sigma_arr_dict[xim_str], 4, 4, '-'))
print(lsq_pm44)
m_44test = Minuit(lsq_pm44, Om_ch2=0.105, AS=4.05)

m_44test.limits['Om_ch2'] = (0.03, 0.7)
m_44test.limits['AS'] = (1.5, 6)

m_44test.migrad(ncall=25)<jupyter_output>/home/yve/anaconda3/lib/python3.7/site-packages/pyccl/core.py:678: CCLWarning: CCL does not properly compute the linear growth rate in cosmological models with massive neutrinos!
  category=CCLWarning)<jupyter_text>Fitting to all 20 correlation function datasetsThis fits for the cosmological parameters using all 20 correlation functions with the random Gaussian noise datasets.**Warning:** This is Very Slow (i.e. each cell below takes ~15 min to >1 hr to run).<jupyter_code>#Joint-fitting all auto- and cross-correlation functions with random Gaussian noise for each corr. fctn
lsq_tot = None
for i in range(4):
    for j in range(i,4,1):
        pair_str = str(i+1)+str(j+1)
        xip_str = 'xip_'+pair_str
        xim_str = 'xim_'+pair_str        
        lsq_ij = (LeastSquares(xi_pm_ij, theta_sample, random_noisy_data_dict[xip_str], 
                               sigma_arr_dict[xip_str], i+1, j+1, '+')
                 +LeastSquares(xi_pm_ij, theta_sample, random_noisy_data_dict[xim_str],
                               sigma_arr_dict[xim_str], i+1, j+1, '-'))
        if lsq_tot is not None:
            lsq_tot += lsq_ij
        else:
            lsq_tot = lsq_ij
mtot = Minuit(lsq_tot, Om_ch2=0.105, AS=4.05)

mtot.limits['Om_ch2'] = (0.03, 0.7)
mtot.limits['AS'] = (1.5, 6)

mtot.migrad()
mtot.draw_profile('Om_ch2', bound=2, size=20)
mtot.draw_profile('AS', bound=2, size=20)
mtot.draw_contour('Om_ch2', 'AS', size=10)<jupyter_output>/home/yve/anaconda3/lib/python3.7/site-packages/pyccl/core.py:678: CCLWarning: CCL does not properly compute the linear growth rate in cosmological models with massive neutrinos!
  category=CCLWarning)
{'fraction_non_alphanumeric': 0.13838262577758376, 'fraction_numerical': 0.1364182036450944, 'mean_word_length': 3.431726583911011, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 114, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '403420', 'license': 'MIT', 'n_tokens_mistral': 13039, 'n_tokens_neox': 10726, 'n_words': 2088, 'path': 'pyccl_shear_fitting.ipynb', 'repo_name': 'cahebert/psf-to-cosmo', 'size': '27504'}