diff --git a/fgspectra/cross.py b/fgspectra/cross.py index 0acd332..4262449 100644 --- a/fgspectra/cross.py +++ b/fgspectra/cross.py @@ -8,7 +8,7 @@ import numpy as np from . import frequency as fgf from . import power as fgp -from .model import Model +from .model import Model, _apply class Sum(Model): @@ -176,28 +176,21 @@ def diff(self, sed_kwargs={}, cl_kwargs={}): diff : dict dict with same keys as the parameters passed to diff which stores derivatives with respect to parameters """ - sed_diff = self._sed.diff(**sed_kwargs) - sed = self._sed(**sed_kwargs) #shape of sed is ``(..., freq)`` - cl_diff = self._cl.diff(**cl_kwargs) - cl = self._cl(**cl_kwargs) #shape of cls is ``(..., ell)`` - tot_diff_sed = {} - tot_diff_cl = {} - for param in sed_diff.keys(): - if sed_diff[param] is None: - tot_diff_sed[param] = None - else: - diff = sed_diff[param] #shape of diff is ``(param,...,freq)`` - tot_diff_sed[param] = np.einsum('...i,p...j,...l->p...ijl', sed, diff, cl) + \ - np.einsum('p...i,...j,...l->p...ijl', diff, sed, cl) - - for param in cl_diff.keys(): - if cl_diff[param] is None: - tot_diff_cl[param] = None - else : - diff = cl_diff[param] # shape of diff is ``(param,...,ell)`` - tot_diff_cl[param] = np.einsum('...i,...j,p...l->p...ijl',sed, sed, diff) - - return {'sed_kwargs':tot_diff_sed, 'cl_kwargs':tot_diff_cl} + sed = self._sed(**sed_kwargs) # shape of sed is (..., freq) + cl = self._cl(**cl_kwargs) # shape of cls is (..., ell) + + def diff_in_sed(diff): + res = np.einsum('...i,p...j->p...ij', sed, diff) + res += np.einsum('p...ij->p...ji', res) + return res[..., None] * cl + + def diff_in_cl(diff): + return np.einsum('...i,...j,p...l->p...ijl',sed, sed, diff) + + return { + 'sed_kwargs': _apply(diff_in_sed, self._sed.diff(**sed_kwargs)), + 'cl_kwargs': _apply(diff_in_cl, self._cl.diff(**cl_kwargs)) + } class CorrelatedFactorizedCrossSpectrum(FactorizedCrossSpectrum): @@ -249,9 +242,41 @@ def eval(self, sed_kwargs={}, cl_kwargs={}): Cross-spectrum. The shape is ``(..., freq, freq, ell)``. """ - f_nu = self._sed(**sed_kwargs) + sed = self._sed(**sed_kwargs) # shape of sed is (comp, ..., freq) return np.einsum('k...i,n...j,...knl->...ijl', - f_nu, f_nu, self._cl(**cl_kwargs)) + sed, sed, self._cl(**cl_kwargs)) + + def diff(self, sed_kwargs={}, cl_kwargs={}): + """Compute the derivative of the model with respect to every + parameters. + + Parameters + ---------- + sed_kwargs : dict + Arguments for which the `sed` is evaluated. + cl_kwargs : dict + Arguments for which the `cl` is evaluated. + + Returns + ------- + diff : dict + dict with same keys as the parameters passed to diff which stores derivatives with respect to parameters + """ + sed = self._sed(**sed_kwargs) # shape of sed is (comp, ..., freq) + cl = self._cl(**cl_kwargs) # shape of cls is (..., comp, comp, ell) + + def diff_in_sed(diff): + res = np.einsum('pk...i,n...j,...knl->p...ijl', diff, sed, cl) + res += np.einsum('pk...i,n...j,...nkl->p...ijl', diff, sed, cl) + return res + + def diff_in_cl(diff): + return np.einsum('k...i,n...j,p...nkl->p...ijl', sed, sed, diff) + + return { + 'sed_kwargs': _apply(diff_in_sed, self._sed.diff(**sed_kwargs)), + 'cl_kwargs': _apply(diff_in_cl, self._cl.diff(**cl_kwargs)) + } class PowerLaw(FactorizedCrossSpectrum): @@ -336,7 +361,7 @@ def eval(self, nu=None, ell=None, nwhite=None): print('Got {:d} white noise levels, expected {:d}'.format( len(nwhite), n_freqs)) res = np.broadcast_to(np.diag(nwhite**2), (n_ell, n_freqs, n_freqs)) - return np.transpose(res, (1, 2, 0)) + return np.einsum('lij,l->ijl', res, ell*(ell+1)/2./np.pi) def diff(self, nu=None, ell=None, nwhite=None): """ Evaluation of the derivative of the model diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index 1a94602..b6457f3 100644 --- a/fgspectra/frequency.py +++ b/fgspectra/frequency.py @@ -11,7 +11,7 @@ import numpy as np from scipy import constants -from .model import Model +from .model import Model, _apply from functools import wraps @@ -69,9 +69,10 @@ def eval(self, nu=None, beta=None, nu_0=None): nu_0 = np.array(nu_0)[..., np.newaxis] return (nu / nu_0)**beta * (_rj2cmb(nu) / _rj2cmb(nu_0)) - def diff(self, nu=None, beta=None, nu_0=None): + def diff(self, **kwargs): """ Evaluation of the first derivative of the SED + Parameters ---------- nu: float or array @@ -89,11 +90,25 @@ def diff(self, nu=None, beta=None, nu_0=None): Each key of the dict corresponds to a parameter of the model. """ - (nu, beta, nu_0) = self._replace_none_args((nu, beta, nu_0)) - beta = np.array(beta)[..., np.newaxis] - nu_0 = np.array(nu_0)[..., np.newaxis] - sed_diff_beta = beta*(nu / nu_0)**(beta-1.) * (_rj2cmb(nu) / _rj2cmb(nu_0)) - return {'nu':None, 'beta':sed_diff_beta[np.newaxis, ...], 'nu_0':None} + if 'nu' in kwargs or 'nu_0' in kwargs: + raise NotImplementedError( + 'Derivatives with respect to nu and nu_0 are not implemented') + + defaults = self.defaults + if defaults['beta'] is not None: + return {} + + beta = np.array(kwargs['beta']) + nu = defaults['nu'] + nu_0 = defaults['nu_0'] + res = np.zeros((beta.size, beta.size, nu.size)) + + np.einsum('bbf->bf', res)[:] = ( + np.log(nu / nu_0) + * (nu / nu_0) ** (beta.reshape(-1, 1)) + * (_rj2cmb(nu) / _rj2cmb(nu_0)) + ) + return {'beta': res.reshape((beta.size,)+beta.shape+nu.shape)} class Synchrotron(PowerLaw): @@ -131,7 +146,7 @@ def eval(self, nu=None, sed=None): return None return np.array(sed) - def diff(self, nu=None, sed=None): + def diff(self, **kwargs): """ Evaluation of the first derivative of the SED. Parameters @@ -146,17 +161,15 @@ def diff(self, nu=None, sed=None): sed_diff: dict Each key of the dict corresponds to a parameter of the model. """ - (nu, sed) = self._replace_none_args((nu, sed)) - if type(sed) in (int, float): - sed = [sed] - if type(nu) in (int, float): - nu = [nu] - try: - assert len(nu) == len(sed) - except AssertionError: - print('Size of SED must match number of frequencies') - return None - return {'nu':None, 'sed':np.eye(len(sed))} + if 'nu' in kwargs: + raise NotImplementedError( + 'Derivative with respect to nu does not make sense here') + + defaults = self.defaults + if defaults['sed'] is not None: + return {} + sed = np.array(kwargs['sed']) + return {'sed': np.eye(sed.size)} class ModifiedBlackBody(Model): @@ -191,9 +204,75 @@ def eval(self, nu=None, nu_0=None, temp=None, beta=None): temp = np.array(temp)[..., np.newaxis] x = 1e+9 * constants.h * nu / (constants.k * temp) x_0 = 1e+9 * constants.h * nu_0 / (constants.k * temp) - res = (nu / nu_0)**(beta + 1.0) * np.expm1(x_0) / np.expm1(x) + res = (nu / nu_0) ** (beta + 1.0) * np.expm1(x_0) / np.expm1(x) return res * (_rj2cmb(nu) / _rj2cmb(nu_0)) + def diff(self, **kwargs): + """ Evaluation of the first derivative of the SED. + + Parameters + ---------- + nu: float or array + Frequency in GHz. + beta: float or array + Spectral index. + temp: float or array + Dust temperature. + nu_0: float + Reference frequency in Hz. + + Returns + ------- + sed_diff: dict + Each key of the dict corresponds to a parameter of the model. + """ + if 'nu' in kwargs: + raise NotImplementedError( + 'Derivative with respect to nu does not make sense here') + if 'nu_0' in kwargs: + raise NotImplementedError( + 'Derivative with respect to nu does not make sense here') + + defaults = self.defaults + nu = defaults['nu'] + nu_0 = defaults['nu_0'] + res = {} + + beta = defaults['beta'] + if beta is None: + beta = np.asarray(kwargs['beta']) + if defaults['temp'] is None: + temp = np.asarray(kwargs['temp']) + else: + temp = defaults['temp'] + x = 1e+9 * constants.h * nu / (constants.k * temp) + x_0 = 1e+9 * constants.h * nu_0 / (constants.k * temp) + res_beta = np.zeros((beta.size, beta.size, nu.size)) + np.einsum('aai->ai', res_beta)[:] = np.log(nu / nu_0) * ( + nu / nu_0) ** (beta + 1.0) * np.expm1( + x_0) / np.expm1(x) * (_rj2cmb(nu) / _rj2cmb(nu_0)) + res['beta'] = res_beta.reshape( + (beta.size,) + beta.shape + nu.shape) + + temp = defaults['temp'] + if temp is None: + temp = np.asarray(kwargs['temp']) + if defaults['beta'] is None: + beta = np.asarray(kwargs['beta']) + else: + beta = defaults['beta'] + x = 1e+9 * constants.h * nu / (constants.k * temp) + x_0 = 1e+9 * constants.h * nu_0 / (constants.k * temp) + res_temp = np.zeros((temp.size, temp.size, nu.size)) + np.einsum('aai->ai', res_temp)[:] = (nu / nu_0) ** ( + beta + 1.0) * (x * (np.expm1(x_0) * np.exp(x)) / ( + temp * np.expm1(x) ** 2) - x_0 * np.exp(x_0) / ( + temp * np.expm1(x))) * (_rj2cmb( + nu) / _rj2cmb(nu_0)) + res['temp'] = res_temp.reshape( + (temp.size,) + temp.shape + nu.shape) + return res + class CIB(ModifiedBlackBody): """ Alias of :class:`ModifiedBlackBOdy` @@ -225,6 +304,9 @@ def eval(self, nu=None, nu_0=None): """ return ThermalSZ.f(nu) / ThermalSZ.f(nu_0) + def diff(self, **kwargs): + return {} + class FreeFree(Model): r""" Free-free @@ -263,14 +345,65 @@ def eval(self, nu=None, EM=None, Te=None): - Free-free emission in temperature. """ - EM = np.array(EM)[..., np.newaxis] - Te = np.array(Te)[..., np.newaxis] - Teff = (Te / 1.e3)**(1.5) + EM = np.array(EM)[..., np.newaxis] + Te = np.array(Te)[..., np.newaxis] + Teff = (Te / 1.e3) ** (1.5) nuff = 255.33e9 * Teff - gff = 1. + np.log(1. + (nuff / nu)**(np.sqrt(3) / np.pi)) + gff = 1. + np.log(1. + (nuff / nu) ** (np.sqrt(3) / np.pi)) print("warning: I need to check the units on this") return EM * gff + def diff(self, **kwargs): + """ Evaluation of the first derivative of the SED. + Parameters + ---------- + nu: float or array + Frequency in GHz. + EM: float or array + Emission measure in cm^-6 pc (usually around 300). If array, the shape is ``(...)``. + Te: float or array + Electron temperature (typically around 7000). If array, the shape is ``(...)``. + + Returns + ------- + sed_diff: dict + Each key of the dict corresponds to a parameter of the model. + """ + if 'nu' in kwargs: + raise NotImplementedError( + 'Derivative with respect to nu does not make sense here') + + defaults = self.defaults + nu = defaults['nu'] + res = {} + + Te = defaults['Te'] + if Te is None: + Te = kwargs['Te'] + + EM = defaults['EM'] + if EM is None: + EM = np.asarray(kwargs['EM']) + res['EM'] = self.eval(Te=Te, EM=1.0)[None] + + if defaults['Te'] is None: + Te = np.asarray(kwargs['Te']) + Teff = (Te / 1.e3) ** 1.5 + nuff = 255.33e9 * Teff + nuff_diff = 255.33e9 * 1.5 * (Te / 1.e3) ** 1.5 / Te + res_Te = np.zeros((Te.size, Te.size, nu.size)) + + np.einsum('aai->ai', res_Te)[:] = EM * nuff_diff * np.sqrt(3) * ( + nuff / nu) ** (np.sqrt(3) / np.pi) / ( + np.pi * nuff * ( + nuff / nu) ** ( + np.sqrt( + 3) / np.pi) + np.pi * nuff) + + res['Te'] = res_Te.reshape((Te.size,) + Te.shape + nu.shape) + + return res + class ConstantSED(Model): """Frequency-independent component.""" @@ -338,3 +471,26 @@ def eval(self, kwseq=None): for i in range(len(seds)): res[i] = seds[i] return res + + def diff(self, kwseq=None): + """Compute the SED with the given frequency and parameters. + + *kwseq + The length of ``kwseq`` has to be equal to the number of SEDs + joined. ``kwseq[i]`` is a dictionary containing the keyword + arguments of the ``i``-th SED. + """ + if kwseq: + diffs = [s.diff(**kw) for kw, s in zip(kwseq, self._seds)] + else: + diffs = [s.diff() for s in self._seds] + n_comp = len(diffs) + for i in range(n_comp): + def expand(sed): + factor = np.zeros(n_comp) + factor[i] = 1. + return sed[..., None, :] * factor[..., None] + + diffs[i] = _apply(expand, diffs[i]) + + return {'kwseq':diffs} diff --git a/fgspectra/model.py b/fgspectra/model.py index f3f5a82..f46194f 100644 --- a/fgspectra/model.py +++ b/fgspectra/model.py @@ -5,6 +5,13 @@ import numpy as np from copy import deepcopy +# Conventions +# * all the Models have an eval method that evaluates the model +# * a parameters that is defaulted to None means that it is a free parameter +# * The output fo the diff method ahs the same shape as that of eval but with +# one and only one extra dimension (the idea is that you use diff with +# minimizers, which work on flat array + class Model(ABC): """ Abstract class for model definition @@ -235,3 +242,13 @@ def diff_kwargs2array(self, diff_kwargs): val = val[p] res_list.append(np.array(val)) return np.concatenate(res_list) + + +def _apply(func, list_tuple_dict_or_value, **kwargs): + if isinstance(list_tuple_dict_or_value, (list, tuple)): + res = [_apply(func, v, **kwargs) for v in list_tuple_dict_or_value] + elif isinstance(list_tuple_dict_or_value, dict): + res = {k: _apply(func, v, **kwargs) for k, v in list_tuple_dict_or_value.items()} + else: + return func(list_tuple_dict_or_value, **kwargs) + return res diff --git a/fgspectra/power.py b/fgspectra/power.py index e2410e0..37558ad 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -11,7 +11,7 @@ import os import pkg_resources import numpy as np -from .model import Model +from .model import Model, _apply def _get_power_file(model): @@ -81,23 +81,58 @@ def __init__(self, filenames, **kwargs): def eval(self, ell=None, ell_0=None, amp=1.0): """Compute the power spectrum with the given ell and parameters.""" - return amp * self._cl[..., ell] / self._cl[..., ell_0, np.newaxis] + return amp * self._cl[..., np.rint(ell).astype('int')] / \ + self._cl[..., np.rint(ell_0).astype('int'), np.newaxis] + + def diff(self, **kwargs): + """ + Parameters + ---------- + ell: float or array + Multipole + ell_0: float + Reference ells + amp: float or array + Amplitude, + + Returns + ------- + cl_diff: dict + Each key of the dict corresponds to a parameter of the model. + """ + if 'ell' in kwargs or 'ell_0' in kwargs: + raise NotImplementedError( + 'Derivatives with respect to ell and ell_0 are not implemented') + + defaults = self.defaults + + if defaults['amp'] is not None: + return {} + + amp = np.asarray(kwargs['amp']) + ell = defaults['ell'] + ell_0 = defaults['ell_0'] + res = np.zeros((amp.size, amp.size, ell.size)) + + np.einsum('aal->al', res)[:] = self.eval(ell=ell, ell_0=ell_0, amp=1) + res_amp = res.reshape((amp.size,) + amp.shape + ell.shape) + return {'amp': res_amp} class tSZ_150_bat(PowerSpectrumFromFile): """PowerSpectrum for Thermal Sunyaev-Zel'dovich (Dunkley et al. 2013).""" - def __init__(self): + def __init__(self, **kwargs): """Intialize object with parameters.""" - super().__init__(_get_power_file('tsz_150_bat')) + super().__init__(_get_power_file('tsz_150_bat'), **kwargs) class kSZ_bat(PowerSpectrumFromFile): """PowerSpectrum for Kinematic Sunyaev-Zel'dovich (Dunkley et al. 2013).""" - def __init__(self): + def __init__(self, **kwargs): """Intialize object with parameters.""" - super().__init__(_get_power_file('ksz_bat')) + super().__init__(_get_power_file('ksz_bat'), **kwargs) class PowerLaw(Model): @@ -129,7 +164,7 @@ def eval(self, ell=None, alpha=None, ell_0=None, amp=1.0): amp = np.array(amp)[..., np.newaxis] return amp * (ell / ell_0)**alpha - def diff(self, ell=None, alpha=None, ell_0=None, amp=1.0): + def diff(self, **kwargs): """ Parameters ---------- @@ -147,12 +182,42 @@ def diff(self, ell=None, alpha=None, ell_0=None, amp=1.0): cl_diff: dict Each key of the dict corresponds to a parameter of the model. """ - (ell, alpha, ell_0, amp) = self._replace_none_args((ell, alpha, ell_0, amp)) - alpha = np.array(alpha)[..., np.newaxis] - deriv_alpha = alpha*amp*(ell/ell_0)**(alpha-1.) - deriv_amp = (ell/ell_0)**alpha - return {'ell':None, 'alpha':deriv_alpha[np.newaxis, ...], - 'ell_0': None, 'amp':deriv_amp[np.newaxis, ...]} + if 'ell' in kwargs or 'ell_0' in kwargs: + raise NotImplementedError( + 'Derivatives with respect to ell and ell_0 are not implemented') + + defaults = self.defaults + res = {} + + alpha = defaults['alpha'] + if alpha is None: + alpha = kwargs['alpha'] + + amp = defaults['amp'] + if amp is None: + amp = np.asarray(kwargs['amp']) + res['amp'] = self.eval(alpha=alpha, amp=1.0)[None] + elif type(amp) in [int,float]: + amp = np.array([amp]) + + if amp.size != 1 or amp.ndim > 1: + raise NotImplementedError('amp has to be a scalar for now') + + if defaults['alpha'] is None: + alpha = np.asarray(kwargs['alpha']) + ell = defaults['ell'] + ell_0 = defaults['ell_0'] + res_alpha = np.zeros((alpha.size, alpha.size, ell.size)) + + np.einsum('aal->al', res_alpha)[:] = ( + amp * np.log(ell / ell_0) + * (ell / ell_0)**(alpha.reshape(-1, 1)) + ) + res['alpha'] = res_alpha.reshape( + (alpha.size,) + alpha.shape + ell.shape) + + return res + class FreeCls(Model): """ @@ -185,7 +250,7 @@ def eval(self, ell=None, cls=None): res = cls return res - def diff(self, ell=None, cls=None): + def diff(self, **kwargs): """ Evaluation of the first derivative of the power spectrum. Parameters @@ -200,16 +265,20 @@ def diff(self, ell=None, cls=None): cl_diff: dict Each key of the dict corresponds to a parameter of the model. """ - (ell, cls) = self._replace_none_args((ell, cls)) - if type(ell) in (float, int): - ell = [ell] - if type(cls) in (float, int): - cls = [cls] - try: - assert len(ell) == len(cls) - except AssertionError: - print('Cls must have same size as ells') - return {'ell':None, 'cls':np.eye(len(ell))} + if 'ell' in kwargs: + raise NotImplementedError( + 'Derivative with respect to ell does not make sense here') + + defaults = self.defaults + if defaults['cls'] is not None: + return {} + cls = np.asarray(kwargs['cls']) + + # Convention: avoid using eye (wasteful for large number of multipoles) + # Reverted to old method for now + #return {'cls': np.ones((cls.size, 1))} + return {'cls': np.eye(cls.size)} + class CorrelatedPowerLaws(PowerLaw): @@ -283,11 +352,11 @@ def set_defaults(self, **kwargs): @property def defaults(self): - return {'kwseq': [ps.defaults for ps in self.power_spectra]} + return {'kwseq': [ps.defaults for ps in self._power_spectra]} def _get_repr(self): return {type(self).__name__: - [ps._get_repr() for ps in self.power_spectra]} + [ps._get_repr() for ps in self._power_spectra]} def eval(self, kwseq=None): @@ -379,8 +448,11 @@ def eval(self, kwseq=None): The length of `argss` has to be equal to the number of SEDs joined. ``kwseq[i]`` is the argument list of the ``i``-th SED. """ - spectra = np.array( - [ps(**kwargs) for ps, kwargs in zip(self._power_spectra, kwseq)]) + if kwseq: + spectra = np.array([ps(**kwargs) for ps, kwargs + in zip(self._power_spectra, kwseq)]) + else: + spectra = np.array([ps() for ps in self._power_spectra]) res = np.empty( # Shape is (..., comp, comp, ell) spectra.shape[1:-1] + (self.n_comp, self.n_comp) + spectra.shape[-1:]) @@ -396,6 +468,29 @@ def eval(self, kwseq=None): assert i_corr == len(spectra) return res + def diff(self, kwseq=None): + """"Compute the first derivative of the cls.""" + if kwseq is None: + kwseq = self.defaults['kwseq'] + def diff_in_cls(diff, i, j): + shape = diff.shape + res = np.zeros((shape[0], self.n_comp, self.n_comp, shape[-1])) + res[..., i, j, :] = diff + res[..., j, i, :] = res[..., i, j, :] + return res + diffs = [] + i_corr = 0 + for k_off_diag in range(0, self.n_comp): + for el_off_diag in range(self.n_comp - k_off_diag): + i = el_off_diag + j = el_off_diag + k_off_diag + diffs.append(_apply(diff_in_cls, self._power_spectra[i_corr].diff(**kwseq[i_corr]), i=i, j=j)) + i_corr += 1 + return {'kwseq': diffs} + + + + class SZxCIB_Reichardt2012(PowerSpectraAndCorrelation): """PowerSpectrum for SZxCIB (Dunkley et al. 2013).""" diff --git a/notebooks/ACT_Models.ipynb b/notebooks/ACT_Models.ipynb index ee53435..50d5248 100644 --- a/notebooks/ACT_Models.ipynb +++ b/notebooks/ACT_Models.ipynb @@ -68,9 +68,50 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "FactorizedCrossSpectrum:\n", + "- ConstantSED (SED):\n", + " amp: 1.0\n", + " nu: null\n", + "- kSZ_bat (Cl):\n", + " amp: 1.0\n", + " ell: null\n", + " ell_0: null\n", + "\n", + "CorrelatedFactorizedCrossSpectrum:\n", + "- Join (SED):\n", + " - ThermalSZ:\n", + " nu: null\n", + " nu_0: null\n", + " - CIB:\n", + " beta: null\n", + " nu: null\n", + " nu_0: null\n", + " temp: null\n", + "- SZxCIB_Addison2012 (Cl):\n", + " - PowerSpectrumFromFile:\n", + " amp: 1.0\n", + " ell: null\n", + " ell_0: null\n", + " - PowerLaw:\n", + " alpha: null\n", + " amp: 1.0\n", + " ell: null\n", + " ell_0: null\n", + " - PowerSpectrumFromFile:\n", + " amp: 1.0\n", + " ell: null\n", + " ell_0: null\n", + "\n" + ] + } + ], "source": [ "# define the models from fgspectra\n", "ksz = fgc.FactorizedCrossSpectrum(fgf.ConstantSED(), fgp.kSZ_bat())\n", @@ -86,7 +127,10 @@ "\n", "# for testing purposes we'll also compute the tSZ and clustered CIB alone\n", "tsz = fgc.FactorizedCrossSpectrum(fgf.ThermalSZ(), fgp.tSZ_150_bat())\n", - "cibc = fgc.FactorizedCrossSpectrum(fgf.CIB(), fgp.PowerLaw())" + "cibc = fgc.FactorizedCrossSpectrum(fgf.CIB(), fgp.PowerLaw())\n", + "\n", + "print(ksz)\n", + "print(tSZ_and_CIB)" ] }, { @@ -100,7 +144,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -110,7 +154,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 16, "metadata": {}, "outputs": [ { @@ -120,7 +164,7 @@ " 3.49300350e+00, 3.49533489e+00, 3.49766706e+00]]])" ] }, - "execution_count": 8, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -149,7 +193,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ @@ -197,7 +241,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ @@ -206,22 +250,38 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 10, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0.01754513 0.02625792 0.03452104 ... 9.79613729 9.79752592 9.79892278]]]\n" + ] + } + ], "source": [ "# remove the last tSZ, CIBC since they are redundant\n", - "specs = specs[:-2]" + "specs = specs[2]\n", + "print(specs)" ] }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 19, "metadata": {}, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "No handles with labels found to put in legend.\n" + ] + }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEmCAYAAADP6P+fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmcW3d97//XV/s2MxrN2B7v47Ed29ljO8EOEJJgh9zShHAZJwQItNDYpbRwC218Sfsr0NKmzr2Uy4X7KDYlpAFCgh22piHEplmakDRekthZ7Dge7/bYs2kW7Tr6/v44kkaaXTOakTT+PB+P8zhHXx1JX0m23vP9nu/5HqW1RgghhCglS6krIIQQQkgYCSGEKDkJIyGEECUnYSSEEKLkJIyEEEKUnISREAVSSq1USjWVuh4TMR3eg5hebKWugBDlRCnlB7ZrrdePsNsmYKNSqgVoSZcFAD+wWWu9Y4jnbQI2A0fSRYvT+waLVvn+16r49yAuPBJGQqQppdYBW8e4+y5gNdCE+WO+C9iitW4ZuKNSaiXwPa31qpyyJmCvUmpVMX/Mp8N7EBcmCSNxwVNKbQFWAtuBfentkQRHaXUMtB2zJZKltW5RSm0FvgdsKOC5hjQd3oO4sEkYiQue1npzZlsptWqkfQuVbqk0aa13DXH3DmCLUso/0ZbFdHgP4sImAxiEmFwb6D8mkyenO+z2qavOuEyH9yDKnISREJNrNcP8kKcFgaK2ZCbBdHgPosxJN50Q46SU2og5+iyIObJs6xAH/zODA4bTmd4HpdR2YF36OQHuH9D9tpf+Y0EtWuvF5fYehBgvCSMhxkEptVFrvS3nth9zZNkGrfW+Ap8uAKC13pB+rp2YoXTfgP3uA7YA64ca8VaoyXgPQoyXdNMJUbituT/iAOmD91sxR53l8jO6vH3So9xaMIMn13qt9eJiBBGT/B6EKJSEkRAFGiEMdgFNSqnmIrzMBsyTUpsh25021vOHRjVF70GIMZMwEqJ4Mj/wVxf4uEFDotPdZPcD30sPrWYcXWfjUbT3IEQh5JiREEWitQ4qpSD/hNMgo3dhdQ7zfJvTQbRda11bnFqOrNjvQYixkpaREAVQSm1PDzAYSW4rYQ8jH9wfbaTaLsCf7qYrihK8ByFGJWEkRGFWMvrIsd052/sYvVUxZDCkj9vsxJyGZ2sRZ9mesvcgxFhJGAlRmB25k4Xmyjnonzttzk6GOQcnPfnowP0z9/lJT8GTHvW2g8Gj3MZrSt6DEIWQMBKiMDszAwqGsAnzhz470CA9n9u+YUan3ZHef6iD/1/WWt+fc/tuzFFuA4d7j8dUvQchxkzCSIh8AUbowkr/MK8f+MOcDokAZmgMdDfpyURz9m8Cmgfur5Typ2di6BjwukHMk17vGcOw65K+ByHGQ2mtS10HIUoqPThgA2ZXVKY7Koh5UH6P1nrTMI9ZRf8F6XYOaMkM3H/UC9PlzLwAQ0z3o5Q6klO/fcB9mYvglct7EGK8JIyEEEKUnHTTCSGEKDkJIyGEECUnYSSEEKLkJIyEEEKUnISREEKIkpOJUtOUUjKsUAghxklrrSbyeGkZCSGEKDlpGQ0g510JIcTYpS85MmHSMhJCCFFyEkZCCCFKTsJICCFEyUkYCSGEKDkJIyGEECUno+km6MEHH2Tfvn2EQiFCoRB/8Rd/werVq0tdLSGEqCgSRhP0+OOP89hjj2Vvf+QjH5EwEkKIAkk33QR5vd682319fSWqiRBCVC4JowkaGEahUKhENRFCiMolYTRBEkZCCDFxEkYT5PP58m5LGAkhROEkjCZIjhkJIcTESRhNkHTTCSHExEkYTZB00wkhxMRJGE1QVVVV3u1gMFiimgghROWSMJqgWbNm5d0+d+5ciWoihBCVS8JoggaGUWtra4lqIoQQlUvJlU1NSikNhV/pNRQK5R03cjgcRKPRol39UAghylnmt05rPaEfPWkZTZDX680Lo3g8TkdHRwlrJIQQlUfCqAgaGxvzbh8+fLg0FRFCiAolYVQEy5Yty7t96NChEtVECCEqk4RRESxfvjzv9ltvvVWimgghRGWSMCqCiy++OO/2yy+/XKKaCCFEZZLRdGnjHU0HcPToUZqamrK33W433d3d2O324lVQCCHKkIymKyONjY3Mnj07ezsSifDiiy+WsEZCCFFZ5LLjRaCU4sYbb+THP/5xtuxnP/sZ1113XQlrJYS4kGmtiUajhEIhwuEwoVAob3vgOhQKsXTpUu68886S1Fe66dIm0k0HZvh85CMfyd6eM2cOx48fx2aTvBdCDJZKpYhEIgWFxVjKcteF/p7deuut/PKXvyzoMcXqppNfyiK5+eab8Xg8hMNhAM6cOcNjjz3GHXfcUeKaCSHGwzCMogXDUGWRSKTUb3GQUl51oGzDSCm1U2u9vtT1GCuPx8MnPvEJtm3bli277777aG5uxmq1lrBmQkxvhmEQCoXo6+sbtAxXPtw+uWERi8VK/damXCnDqOy66ZRS64AmYOtEm30Fvu6EuukADh48yIoVK/LKvvvd77Jp06aJVU6IaUBrTSwWGzUcCg2RcmxhlAun04nH48Hr9eL1erPbw5UtWrSIz3zmMwW9RrG66SYljJRSzUBQa71riPtWAquBFszQaRlmvy6tdW3RKzeMYoQRwJ133skjjzySve3z+di3bx9Lly6dWAWFmGKGYdDb20tvby89PT1DrgsNkFQqVeq3VVbcbveIYVFIkAy8z+12T8kx67INo3TLZjuwYWDIKKUyLZ71OWXbgc1a65YB+1ZkGJ08eZLly5dnjx0BXHLJJTz33HMEAoGJVVKIURiGQV9f34gB0tPTM+J9mXXuv+ELkVIKj8dTlGAYLiwslso/u6bsBjCkg2YzsBfoHGa3TcDWAWVbgS3AhmLVpZTmz5/PN77xDT772c9my9544w1uvvlmnnjiCerr60tYO1GutNaEQiGCwSDd3d2DltzykQLkQr3svc/nyy6ZmfRHWwbulwmJTGi43W65FMwUmqxuuiPApiFaRkeA9bmtIKWUH+gamKqV2jLKPMfHPvaxvO46gKamJn7+859z+eWXT/g1RPnQWhOJRIYMkrGESyZgDMMo9VuZdA6HY0zBMNYA8fl806aFUanKtpsOhg6jTOgAtVrr4ID9NbB4QEhVbBiBOQvDf7v5Zp597rm8crvdzle/+lW+9KUv4XQ6i/JaYuK01vT19dHV1TXkEgwGR7wvHo+X+i1MiqqqKqqqqqiurs5b5y5jDRCv14vD4Sj1WxJFVolh1AQcGarC6SBYP2D/cYeRUuqZcTzsfVC8MALo7e3llt+/hWefe3bQfQsXLuRrX/sad955p/wHLaJoNEpnZ+egpaOjY9gwyQTKdGmZ+Hy+QeExnrXX65UWR5kwUiniKYNEysiuc7dzy5KpFEmdIplKYej+7aTOv89c95fN9/r5SGPhvTaVGEYrgb2jhVF6v3WYx5HuB3YONdpulNd/ZhzVLnoYgXnl1z/57Gf5/gMPDHl/Q0MDGzdu5BOf+ISMuMsRiUQGhclwIZN7u5IPurvdbmpqaqipqcHv92e3B5ZVV1cPGyI+n08CZAha60E/4IlUingqmbce7gc+bhgkdHo9bBCM9nyDy4Z8bp0ibuQ/j2byT8FZN2cpOz9Q+Gko0zaMil6ZMSp2N91AP330p3x20x/T2d017D5XXHEFt912GzfccAPvete7cLlck1KXqZI5ljJaqAws6+joIBqNlrr6BXE4HNmwGCpIxhIw07GFnNIpYoZBzEgSNRJEjWR6MbdjqfTtZCJbHksZOfv2PyY24LHZ5xhUbpbl/sAndWUPKVcoHFYrdmU11xYrDou5tlssOCw27BYLdosVm7Jgs1iwKWt6bclZWwfc7i9fVjODzy6/tvC6VWoYMfwxo2kdRgBtbW38zV/ey7/86EGSRnLEfZ1OJ6tWreKKK67giiuu4PLLL2fx4sXMmDFjykf4GIZBd3d3Niy6urrGvF1JZ7F7PB5qa2vx+/3U1tYOWkYqr4SRVymdIpJMEDYS5jqZIJLZNuJ5ZeFkPH1fctB9ESNhhoZhDB0wObfjqfLp+rQp88c686Oeu535gc/9Yc/9gc/88PcHwChlmedTFhzWkZ5vbK9hLePWbiWGUWYAw8CBCkOWT7WpCKOMd94+zD9s/gqPPvFzwvHCWgBut5vGxkYWLlzIjBkzqKuro66ujkAggMPhwGKx5C2pVIpYLEY8HicWi+VtZ6Y9yZywmLudW1ZJZ7jbbLbs5zHUMlyw+P3+shpQkkgZ9CVi9CXihJJx+pIx83Yynr9OjHDfgHUoWZpBFi6rDafVhstqx2W1pRc7Tos1p8yes1/+voPL7DitQzzW0r+fw2r+kOf+sJf7HwuTKWGk6AwnaA/F6QjF6QzH6Qwn6EivO8NxLprh40vXLy74uSsujHLKN2it9+WUNWF2303ZyLmhTGUYZXScPsd3v/5/+MWv/419J94kVWZTM5WS3W4fMlRGK/P5fCX90TFSKXoSUbriEYKxCMF4hGA8ml7n3+6OR4cNkslqUbitdjw2Ox6bI7ttrh24bbb+7UH32XFbbXhsDlxWc79hQyEnPBwXeAhMhkjCoCMUpyMcp70vTkc6VDJB07/dX94THbknBuD6xXU8/ScXQDddunwrZvBsyylrBu7QWpf0pNdShFFGvDvMW0/+F7989DF2v3OAPcffpLWnfcrrMRkcDsegABlLyHi93pL9iCVSBp2xMB2xEB3RMJ2xsBku6aVrhJDpSRTnWJdVWfDZHfhsziHWTnw2x9DrEfbx2OxYVPl291xotNb0xQwzOMJmkPRvJ4YuDycIxwv/Q8VqUQQ8duo8Duo8duq8DgIeBwGPPVveVOfhpmUzC37ucg+jLuBurfWOAeV+YPuA6YB2YgZXybro0vUoWRhlpOJJOl49RvueFlreeJu3zx/n7fPHeafjFCd7z3O6vZW+cGnOsK+pqcl2c+V2eY22XepjKVpr+pIxzkX6OBfpNdfRXs5FemmPhuiIhdPBkw6fWJjexPiPcykUNQ4XfocLv8Odswy+Xe1wUWV3Dhk0TqtNWhQVRmtNKG7Q1henLRTjfF+ctr4B65zytr440WThAyscVgt1XjNA6r0O6rxmwGS2670OM3Sy23ZqXHYslsn591R2YZQOmi9jTn7ajDkR6i7Modk7cvbLDN3OTJS6r5QDFzLKIYxyRdt76XzlGJ2vHSfW2QeYdetJRuj2Jgk64kQdKXpT0eyAgWQyidaaVCqVXZRSOJ1OHA4HTqczuzgcjrx5sgaenJi7XY5nuPclYpwJ95gBE+3NC5vWbJl5O2IkCnpuq7IQcLoJOD3UOb3UOT0EnJ5hQ6XW6cmWV9md0vqYRsLxZDo84pxPB8j5vhhtoaEDJpIoLFxcNgszfI4BwZLZtg8RLA58zvLq+iy7MKp05RZGGVproud7CL51iuCbpwmfyp/2z+p2UNU0E1/jDLzz6/DMqcViq9zrJ2mt6YiFOR3u5lQoyOlwz4B1N6fC3XQXMPDDZbUxy13FLJfPXKe3Z7h81Lk81OWETp3TS7VDAmU6C8eTtPbGONcby1sPLDvXFyu4SywTLjN9zv6118EMn5OZvpy117zf6yivYBkPCaMiK9cwGijRG6H3yDl6jpyj98g54sH8kzyV1YJnTi3e+XXZcHLW+VBl0rLRWtMeC3G0t5NjfZ3pdRdH+8zbx/u6iI4y7B3AabUxx11NQyZc3L4hA2eWu4oqu7Pi/8OLkcWShhkgvXFae6Mjhk1vbPR/XxlOm4UZXgczq8xQGRwy+WXTIVwKJWFUZJUSRrm01sS7QvQeOUffyQ5CJzqItnUz8GRti92Ku8GPe7YfT4Mf9+xa3A01WJ32SamXkUpxItTFoe42DnW38U5vezZ8jvV1jTrE2O9wM9dTzTyvn7meGuZ5a5jrqckrq3N6Lrj/9BeiSMLgbE+UM91RzvTEONMT5XR35nZ/6HRFxt4V67BaaKh20lDlZJbPSUN1el2Vs13tYqbPQZVz+h2701qTivZi9HVihNJLXwdWbwDfpYVfXFvCqMgqMYyGYkTjhE51EjrRQehUB+GzXSS6hz5PyOH34JpRjau+GteMKnN7RjW2KteY/gPGjCQHu89zoOssB4PnOdRjhs/hnrYRWzd+h5tGXy2LfAEaqwLm2lfLoqoAjb4APnv5nO8jJkfCSHGud3C4nOmJZbdPd0fHHDJWi2KWz8msKgcNVS4zaKrSAZO7Xe2ixjV9AkanDIxQF0Zfh7n0tmP0dZDsax9Uli0PdcIQ/z+9l6xn4T1PFVwHCaMimy5hNJRkOEb4bJDI2SCR1i7CZ4NEz/egjaEPtlqcdlwzqnAGfDgDPhy1Xno9isOWXl6NtbO/u5X9nWc52H1+2GlWZrurWVYzg2U1M7ioegaLsqETwO90T+bbFSUWjic51R3lZFeEk8EoJ7sjnApGzNBJB875vhhj+a9msyjm1LiYU+1iTrXTXNe4mFvjYnaVi9npFk7A45i00WJTRScT6cDowBgQJsnMdqY8EzDhLsb0QQ6gnF6s3gA2Xx0WXwCrN4B74Urqb/ly4c8lYVRc0zmMhqKNFLGuENG2npyll2hbD0Zk+G60BCnO2uOctsc4Y4sTr7Lh9fuora9l9qx6FjfM5qLALKodlT2vnhhaLGlwujvKyWA6aIKRQdud4dFbM0rBLJ+TuZmgyQ2c9PbcGhd1FRwyqXiEZM95jN42kj1tGL3nB6zbzPvTAZOK9IzrdSweP7aqeqy+Oqy+zLoOa1U9tsx2brmvDksRex8kjIrsQgujjGAswp6Ok+xuP8nLbSfZ3XaSUG+IxoSLuQlnenGw0HAzP+mmJj76QAibz4Wjxo2jxoO9xoMjvdir3ebic03a8Soxflpr2kNxjnVGONYV5lhneEDgRDjfN/qUQnarYl6Nm/l+F/P9bub73cyrcTHP704Hj3lcxmYtj0E1Y5WKR7LB0h8y5jqznczZ1rECzwm0WLF6AyMGSX/opMu9tShr0S7YPS4SRkV2oYRRMBbh2dYjPN16hKfPvsP+rrOD9qmyO1ldN49rZizg6vr5rK6fxwJvLUopUokksa4Q8c4Qsc4+4sEQ8e4w8e4I8e4wid4IpEb/DC0OG/YqF3afC1uVGVD2Khf2Knd+udeJqrAfrXKltaatL54OmgjH0+tM8Bzriow6lNlqUcytcTGvpj9ockNnvt/NDG9ltGYyB/KT3a0kg60YPefM7e5Wkt3nSPacywucVLSvoOdXNgfW6pnYqmZgrZqBrXpm/rpqBtbq9NpXh8VdUzajXgshYVRk0zmMjvV28tjx/ew4doCX20/kzYHntNq4MjCHq+vnc029GT4X1dSP+zwbnUqR6I2awZQOKHM7TLwnTKI3SqI3ik6O8fwNZZ5LZfM4sXud2DKLx4nN68LmdZhrT/99Vkdp/1Ispa5wnCMdYVo6MgGTHzijnZTpd9tprHXTGPCwsNbNglp3Xug0VLmwlnnQpGJhM1B6zpEMtma3jdygSW/rQqZvstqxZcJluJDJhEv1TCyuqmkzUGIkEkZFNt3C6EhPO48dP8D2Y6+xp/1UttxusbJmxgJumL2EGxoWs2bGQly2qe0y01qTiiWywZToi5LojZDojZLsi/SX90ZIhmODhqqPRtmteeFldTmweRxY3eZic9uxup3ptQNbutziKP9RVlprzvbEONIR4kh7mHc6QhxpD3GkI8w77aFRR5/lhk1jwE1jrSe7vbDWg99dnt2nWmuMUCfJztMkg2dIBM+Q7DrTHzrdremwOUcq2jvm51VOL7aaBmzVs7D5G/q3axrMwMkJGou7uuz/fZSChFGRTYcweru7jR3H9rPj2H5e6TydLffaHPz+/ItpbryMm+cur6ih09pIkYzESYZi5hKO9W+nbydCUZLh/n3G3OoayKKywZRdexxmmLkdWN12rC47VqcDq8uG1eXA6rRjcZnlFntxTng0UprjXWa4HOkIc6Q9lN1u6QyP2JXmcVhZXOehKeBhUZ2n7MMm21XWdcYMmS4zbJJdOYETNBc9xktgKLszHSZmqJhLzna2fBYWl2+S3+H0J2FUZJUaRm8Fz7Hj2H62H9vPgZzjP1V2J7fMv5jmxsu5ee5y3FPc+ikVrTWphEEyFCUZipMMRTEicZKRRHodH3atExO8bINFYXWmA8tlH3bbkt6OKMWpcIJjvQne6YlxKBjh9Y4Ib3aEiRvD/zus9zpYXOdhcZ2XxfUeltR7ze06D7Oqyme2CW0kzYDpOEmi6xTJrtMkcsIlEzhjPdBv8fix++dgq52DrXYutprZ2PyzB4XNaC0YrTWGYZBMJkmlUiSTSQzDyC7DlWfuG27fQp9rqOcba73G+lqjvX5mDkutNWvWrOGBBx4o+HsuVhhduJ3rFexg8Dw/OfoKO47t583guWx5jcPFh+ZfQnPj5ayfc9GUd7+VA6UUVocNq8OHs8ArZKUSRn9AReMkw/mBZUQSGLEERjSBEY1jxJJ52zpppPcb+0Xs7MDS9PLfMvWwQtSmSFgsaLsV5bDhcNlxexz4fE5cbocZaA4bFqfGmoxiCSaxRsL0Omxma81hw+q0YXGklyLPV6hTBolgK7H2Y8TOHyfWfpxE5wkSnacwgmdIBU9DXxuM4XLfhsVB3FlL3F5D1FZF2FpF2OIlbPHRi5te7aYn5SSS0CR7k8Q74iQSCRKJkyQSLSQSCeJxsyyZTGa3M+WZstwfYDFYKFSaKwJkSMsordxbRm3RPh5peZUfHtnL7vaT2fJah5vbFl5K88LLWTdnKY4SD/O8UCWMFG+39vLWyS6OnO3hxPkezraHCHZH8OgUPsCrdM5aU2NR1NkVfqvCo8CZSmEzUjDMycgToRWklMZQmiQpEtogmTKIp5IkUkniRoKYkSCWjBNJxDCSYVSyD5vRh1334dR9uAnhVSGqrSFqrCFsxEHHQEcBg4F/Fqc0dMUtnI9aaY9ZaY9a6YhZ6IhZ04uF9piVsKFg0KMnl9VqzVtsNhsWiyVvPdb7B+471GPGev94Xn+8z2WxWLBardltu92O1+st+LOUbroiK8cw0lrz3LkWvvPWC/zi+OvZ2Q6q7E6aGy/no4uu5IbZS7BbKneW7kqTNFIc6QjzRmsvr7f28kZrL2+c6+Xttj4SQ3StKQULa90sm+Fj2UwfF9V7WVhtY7YrRY0lQSQSIRwOZy8BHw6HCYfCxEIRYqEIiXCMRCSGEU+SiidJpVtgygCLobFohVUr7Fhx25247E5cNgduuxO3zWmu7U6sk/xvJKUNDCNOIpUgbiSJpQyiKYOETpEkhUGKFJqU0qQUaEt6UYDNgrIqsFpQNivKqrDYbVjs1vRiw+qwYXPasTns2Jx27A4Hdrsdu92Ow+HAZrPhGKUs94e4XLoypwMJoyIrpzCKG0l+eGQv//fN57PnAVmVhQ/MvYhPLl7NrQsuuWCOAZVSRyjOa2d6ePVMN6+e7mH/2R7eOtdHfJiWy0xnilm2OHWEqI534QqdxxI8S1+wg+7u7uySSBR2faWxcrlceDwePB43M6sczPVCgzvFDGecenuSWkuSKuJ4UnGsKLRyQnrp33agLS601QvOADhrUPZqsPlQVg9YXSjsaG2FZIpUwiCVMIadWmqyKJsFi82Kslux2Kz9wZUps9uw2CzmoBKbLe9+s8yKxWZJr639zzfsOr2/1Yoq86HtU03CqMjKIYxiRpIHDr/Mffv/g5OhIAAzXT42LVvDpmVrmeutKVndpivDMOjs6uLVo63sOd7Ba2d7ONQZ51jIQtAYusvTGu7EFjyDrfsM9m5zbetpxWKM7ViRy+WiqqoKr9ebDg9P3vZQt3PLvF4vbqcDt9GLM3wOW98Z6DpJoq2FxPkjxM8fGXV4s8VVhX3GIuz1jTjqG7FnlhnmbYvHX1DrQRspUokkqbiBkV6bt5PZwEolkuiEQSpp9IdYZnuosnT5UI8pJWW1mCFltaLs6XVeuA23HiLkrFYsdjPkhg1Dq7lWNiuW9Gsra3opgxaehFGRlTKMUjrFj4+8wr17n+BUuBuAS/yz+J+X38iGxitwynGggiQSCdra2mhra+PcuXOcP3++f7u9neM9BqeSbtotNfR5ZpLwz0PbB8+lp5IxbMHT2LtOYg+ewh48ha37DJZkDIfDgd/vp6amJm/x+/2DyjO3M2uXa2zz9qUSMRLtx4ife4f4+SPpoElvtx0dcahzscOmnGit8wJrzIGWs49OpsuTqf61kd7PSKUfn0Ib/S2/zH7lJBtK2XAzb2dbcTZLOsDM+3LDLBN2ZrkVZ52PupWLCq+DjKabHva0n+TPXvoFL7UdB+Cy2tn8zZXr+O8LL5OrjQ4QjUY5e/ZsNmCGWzo7O80fLMDw1pMILCRe10gisJBE7Sp0zeAwcMT78Ce6aLCEWeAxWFpjY+kMH4HaBmprVwwKF5drbJfZGI1OpUh2nSJ29hDxs4eItZrr+Lm3SXScGHFGZltNA/aZi3HMXIxj1hIcMxebt2ctweoNVGzYjEYphbLbsNin/udLa402UjkhZpBKpvrXueE1IMTy17mPG2n/zPOm902/tjbyF+Iw0fZi1eJZ4wqjYpEwKpGYkeRrrz7FlgNPk9KaBncV/7jq97hryaoLNoR6e3s5deoUp06d4vTp09ntzNLW1jbi41MOD/G6JhKXvRtmLSVaM4+EbfDlKma5LVw+y8O7GgNcu3gmK+cHmFU1uScCG5Fe4q1vE289ZAZP6yFiZw8Sbz2MjoeHfpCyYK9fmB80M81tx8wmOWGzBJRSKJsVbFaslPa4bTYYc8Mxs22YwZjdzgReOjCz2+lgSyUNnLWFj6QrJummS5vKbrq3u9u4/Zkf8lrnGRSKP7/kvXz1qpuoGqKraDrRWnP+/HmOHTvGsWPHOHr0aHZ98uRJuru7R3y8zWajoaGBhoYGZsycia1+AT2+OZy11NISdXCib/B3N8Pn4Or5/uyyer5/UoMn2ddB7NQbxE6nlzNvEj97iGTwzLCPsVbPxNmwDMfsZTgaluGcnd6ub0TZHJNWVyGKQY4ZFdlUhdGvT73Fnc/+mO54lMXW995RAAAgAElEQVRVdfzrez/Ku2eVrmk8GSKRCG+++SaHDh3KC51jx44RDg/TCsA8sD937lzmzZs3aKmfNZuTCRcvnejmd8c6+d2xLtpD+cdMXDYLVy/ws3ZhbTZ8FtS6J6W7yggF+wPn9OtE09tG97kh91d2J45ZS/vDJh0+zoZlWL3+otdPiKkiYVRkUxFG3zv0Ept+9xgazYcXXMq/XvfRim8NRSIRDhw4wP79+zlw4ACvv/46hw8fxjCG7sH2+/0sWrSIRYsW0djYmF0WLlxIXV1d9h92PJliz8kgTx9p5+l3Ovjdsc5BM07PqnLy7sZa3r0owLsbA1w1twaHrbhdnDoZJ3b6TaInXiV6cj+xU68TO/3GsC0d5fTinHMxrnmX4px7CY45F+Ocsxx73QKUnA8mpiEZwFBhvv3m83z+v34BwN9cuZ6vXLm+Io8NRSIRdu/ezUsvvcRLL73EK6+8Qjye30KxWq2sWLGCFStW0NTUlA2ehQsXUls79Bw9SSPFyyeCPP2OGT7PH+scNCHoJQ1VvHdRgHcvCnBtYy2LAp6itnqMvk6iJ18jevzVdPi8Ruz0m2AMPi9IOdw4Z6/AOfcSnPMuwTnXDB8zdCrvexWi1CSMpsAP39mbDaJvvetDfP7i95a4RmMXCoXYvXs3L774Ii+99BKvvfZa3kmbSilWrFjBVVddxaWXXspll13GihUrcLsHDxwY6ERXmN8cauPJg+fZdbidnmgy7/7lM33csKSOG5bU876mOmYW6ViP1ppk5ykiR/cQPfEK0eOvEjv5mjl6bSClcMxaimvBlTgXXJFt8dhnLJKWjhBFJGE0yf7jzGE+/fyjAPyfa8o/iHp7e3n55Zd56aWXePHFF9m/f39el5vFYuHyyy9nzZo1rF27lmuuuQa/f2zHPCIJg/9s6eDJg208eeg8b53Lv3Lm0npvNnyuX1xHQ3VxujCTPW1Eju4menQPkaO7iRzdPeSxHeVw45x3Ga4FV5rLwitxzbtMRq0JMQUkjCZRa7iHjz77I5I6xZcueR9fuKT8gqi7u5uXX3452/I5cOBA3qzGVquVK6+8krVr17JmzRquueYaqqurx/z853pj/Nsbrfzi9VZ+e7idaM5Jg1VOG+9fWs/Ny2fwgWUzaQx4Jvx+jEgv0WN7iLTsTgfQbhLtxwftZ/HW4m5cjatxZTZ8HA1LpbUjRIlIGE0SrTV/8PyjtEVDvH/2Uu6/+oOlrlLW0aNH+dWvfsWvf/1rXn/99bxBGzabLRs+a9eu5eqrr8bnK6xlcOh8H7983Qygl0505Z23edXcam5ePpObl81kbWMtduv4j69orUl0nCBy+AXCh39H5PALRE/uH3TZAuX04l64ElfT1bgbV+Nuuhr7zMXT9qRQISqRhNEk2X7sNX5z+hABp4eHrvtoyQcr9PT08LOf/YxHH32U/fv3Z8vtdjtXXnkla9as4dprr2XVqlUFTyOvtebA2V4effU0PzvQysHz/d1vTpuFdUvr+dClDdxy8awJdb3pZILoiVcJv/M7Im+/QPid35HsOp2/k9WGa8FK3Iuuxr3oalxNV+Ocs0JaPEKUOQmjSRBNJvjL3Y8DcN+q32OOp3QTnB44cIAf/vCH/PznP8+e4+Pz+fjABz7ALbfcwnve854xDTYYytttfTzyyhkeefV03vGfWred3794Fh+6dBYfWDYTn3N8/8x0MkHk2B7Cbz1D6K2nCR9+YdBsBRZvLZ4l1+JZ+m7cS6/FvehqLM6Jd/cJIaaWhNEk+OGRvZwIBbmsdjafWXrNlL++1ppnn32W73znO7z44ovZ8muvvZa77rqLm266acyTdQ50pjvKj/ae4pFXT/PK6Z5seZ3HTvMVc7j9ijlc1xTANo7uN20kiR5/hdBbTxN662kih58nFc0f5OBouMgMniXX4rno3TgalslQaiGmAQmjIkvpFN9441kAvnz5jVin8IdSa81TTz3FN7/5TQ4cOABAVVUVd9xxB3fddRdLliwZ1/PGkyn+7c1WHnj5JE8ePE8qfQyo2mXjw5c28NGr5vL+pfXjOv4TP3eEvtd/Q9/+JwkfepZUpCfvfsfsZXiX34BnxfV4l1+PrWbWuN6DEKK8SRgV2bOtLRzqbmO+109z4+VT9rp79+7l61//Oi+//DIA9fX13H333Xzyk58saPRbrgNne3jg5RP8aO/p7NQ7dqvitotncdeqedy8fCYue2HHYlLRPkIHn6Fv/5OEXv8N8XPv5N3vmLUEz/Ib8K64Ac/y92GvnTOuugshKkvZhpFSaqfWen2p61Gon7S8AsAnF6+a9MuBG4bBCy+8wEMPPcSvf/1rAAKBAF/4whf4+Mc/Pq5jQQkjxc/2n+U7Lxzj+aOd2fLLZlfx6WsW8PGVc5nhK+zk0/i5I/S+8kt6X/t3Im8/n3cdHou3Ft8l6/Fe9gF8l6zHXje/4DoLISpf2YWRUmod0ASsK3VdCpVIGTx23Owe+2jTlZP2Ou+88w7bt2/nscce4+xZ87LkLpeLu+++m8997nNUVVUV/JytPVG2vXSC7754jLM9McDshvv4yrl85poFrJxXM+ah0Fprosf20rv3F/S+8ktip17vv1Mp3Ivfhfeym/Fd9gHci65GycUDhbjgjflXQCnVDAS11ruGuG8lsBpowQySlqH2G4vM45RSW8bz+FLa3XaSzliYZTUzuLR2dlGfO5lM8pvf/IYf/OAHeYMS5s+fT3NzM5/4xCdoaGgo+HnfOtfL/U8f4cf7TpEwzINBF8/y8afvWcRdq+aNeSScThmE3nqG3r0/o3ffL/OGXFvc1fgu/z2qrroV72U3YfPVFVxPIcT0NqZfmnRr5XvAhiHuawK25HapKaW2K6VatNYtRatpBfiPs+bxj/fPXlq05+zt7eWhhx7iwQcf5MwZc6Zor9fLrbfeSnNzM9dccw2WcQySePlEF/f99h1+8XorABYFt13awJ+9ZxE3LKkbUytIp1JE3nmR7v96hJ7d2/Om2LHVzqVq5YeoWnkb3uXvk+vyCCFGNGIYpYNmM7AX6Bxmt03A1gFlW4EtDBFe09kzrUcAuHH2+Eat5erp6eH73/8+3/ve97IXnWtqauLTn/40GzZsKHhWhIz/bOngK785xNPvdADmSal/ePV8/uL6xSyuH/1kV7MLbh89//UI3f/1KMnOk9n7HLOWUn3N7VSt+jCuxpUyw4EQYsxGDKN0y2YTgFJq8zC7NTM4jPYAOzM3lFIbgcUjvNTO8XbrlQutNXs7TgGwdubCcT9PLBbjgQce4Nvf/nY2hNasWcPnPvc5rr/++nG1ggD2ngzy108e5MmD5qW7q102/uTaRr7w3kVjmhUh2XOe7t/9iOBzDxA7/Ua23BaYT82aj1L9ro/iWniVBJAQYlwmdORYKeXHPEaU12rSWgeVUiilmrTWLVrrbRN5nUpwMhQkGI9Q7/Qy2z2+odRPPfUUX/va1zh27BgAa9eu5Ytf/CLXXnvtuOt16Hwff/3rg+zYbw50qHLa+NL7mvgf1zVR47aP+FhtJOnb/2uCzz1A72uPg2Fe4sFaNSMbQO7Fa+SkUyHEhE10GFMAzPAZ5v4mzEEN097+LvPH/orAnIJbBx0dHdx77708/rg5hdCSJUv46le/yg033DDu+gQjCf72qbf59vNHSaY0LpuFP3vPIjbfuIQ678jHbxJdZ+h6eivBZ7aR7DaPKWGx4rvyFvzXfZqqK35PjgEJIYpqomE0tgvZFCA9Mm8d4E+PqCu4C08p9Uyx6zWag8HzAFziL2yGgN/85jf85V/+JR0dHXi9Xu655x4+9alPYbeP3GoZjpHS/ODlE9z764O09cVRCu5es4Cv3rSMOTXDd8dprYkcfoHOnd+mZ+/Psq0gx+zl+K/7NP5r78LmL3y0nhBCjEXZneChtd4H7APuL3VdCnE81AVAY1VgTPsbhsH999/Pd77zHcCcN+6f/umfmD9//Cd9vn62h8/89DVePmE2VN+zKMD/ve1Srpo3/EStOpmg+78eoePJbxA78ZpZaLFSdXUzgXV/imfZdXIcSAgx6YoSRkop/whddVNOa319oY9RSunR9xresT7zsNlCb+2o+/b19bFp0yaeeeYZrFYr9957Lxs3bhz34IR4MsU//PYw//DbwyQMzdwaF//7lou548rhuwxTsTDB575Px6//d/Zy29aqGdRev5HaG/8Ye2DeuOoihBDjMdEwyhwPCgDZMEoPbMi9f9o73me+/YW+kcOos7OTu+66i1dffZVAIMB3v/td3v3ud4/7dV851c1dP3mFN1p7AfjjtQvZ8vsrqHYN3c1nhLvp3PltOnd+C6O3HTAnI637vXuoWfMxLI7iXOpbCCEKMaEwSo+aa2HwsaMA5mwNF0wYnQqZYTTfO3yXWFdXFxs2bODgwYMsWLCAn/zkJzQ2No7r9bTWfOs/j7L58beIGymW1Hv5l9sv532L64fcPxUL0/nb/0fH4/+IETJbca5Fq6n//S9TtfI2GREnhCipYnTT7cKcCmhfTtnKdPkFwUilCMajAASGubBbJBLhU5/6FAcPHmTp0qU88sgj45q+B+B8b4w/fPRVnnjLHDTxJ9c28r9uWYHHMfjr1Mk4Xc9+n/Zf/R3JoDniz3PRe6m/7St4L75RjgcJIcpCIWEUYOjRc5uB7UDuuUSb0ssFoTsRRaPxO9zYhpipW2vN5z//efbu3cucOXN4+OGHxx1Ee04G+fAPdnOqO0qt284Dd1zBbZcNngdPa03fa//OuYe/SPzcYQBcjauY2fz3eC+9SUJICFFWRpsOyA98GfN8IT+wRSm1HnO49Q7IdtVtVkrdQ/9EqVsupC66jmgIGL5V9MADD/DEE09QVVXFww8/zJw547tGz4/2nuLun75GNJni2sZaHvnEKubXDr5MROzMQVof/nNCB54EzGNCM5v/gapVH5YQEkKUpdGmAwpitnxGlDMc+4LUGQ8DEHAODoY33niDv/u7vwPgG9/4BkuXFj6JaiqlufeJg2x52pyI9e41C/jOhy/DYcs/zpOKR2j7+Vfp+M0/gZHE4qlhxm1fJfD+z6Fs4ztvSQghpkLZnWdUiTpj6TBy5LeMUqkUmzdvJpFIcNddd/HBD36w4OdOGCnu/ulr/OueU9gsim/ddimfvXbhoBZO6OCznH3gj8wrpyqF/31/xMzmv8dWPXP8b0wIIaaIhFER9MTTF6MbMCz64Ycf5pVXXqGhoYG/+qu/Kvh5w/Ekd/xwH4+/eQ6Pw8pjn1rNzcvzw8WI9HL+0Xvoevq7ADjnXsLsT/8LniVrxvluhBBi6kkYFUHUSADgtvZ3hUWjUb75zW8C8Dd/8zcFX301kjC45fu7+Y932gl47DzxR+/iXQvzz2GKHN3D6X++02wNWe3MuOWvqL/lyzJvnBCi4kgYFUE0PY+bK+fy2T/5yU9obW3l4osv5tZbby3o+WJJg//+oBlEDVVOfvvHa7m4oT/MdCpFx5P/xPkd94KRwLngCuZu+hGueZcW5w0JIcQUkzAqgsiAlpFhGPzzP/8zAF/84hcLGsGWNFLc/tBenjzYxgyfY1AQGaEgp7d+nL7XngAgsP7zzLx9i8ycIISoaBJGRTCwZfTss89y+vRpFixYwAc+8IExP4/Wmj/9+ev86o1zBDx2dm3KD6LYmYOc/NaHiLe+jdUbYM7dD1J11S1FfS9CCFEKEkZFkDlm5EoPn3744YcBuPPOOwua/PQbz7Sw9cXjOG0WHv/MNVw+p/8ifb2vPcHpf76TVKQH5/zLmf+FX+CYsaiI70IIIUpHwqgIcltGoVCI3/72tyiluP3228f8HL96vZW/fPxNAH74satY29h/KYrgcz/gzA/uhpRB1eqPMPfuB7G4fEV9D0IIUUoSRkXQH0Z2nn32WeLxOKtWrRrzlD8tHSE++ZNXAPiH31vOhiv6Z2ho//f7Of9T87zj+lvuZcZHvi6zKAghph0JoyLoH9pt46mnngLgpptuGtNjY0mD2x/aS3c0yYcumcX/vHEJYB4/anvs/6P93/4egFkf/xZ1N31+EmovhBClJ2FUBJFk+piR1c7zzz8PwI033jimx977xEH2nuqmMeDmBx+9Mtvqaf/l35lBZLEyd+ND1Kz92ORUXgghyoCEURGsnbmQpE5RE9OcPXuW6upqli9fPurjfne0k28+14LVonj0rlXUesyTVdv/fQttP/8KKAtzN/2ImjUfney3IIQQJSVXVCuCP15+LT9+38extrQCcNVVV406ii6SMPjDR19Fa7jnhsVcs8CcXSH4nw9y/qf/E5Rizh/9QIJICHFBkDAqon37zInLV61aNeq+f7/rMG+3hbh4lo+v3HQRAKE3/8McNQc0fOLb+N/zycmrrBBClBEJoyI6dOgQAJdddtmI+x3tCPO/nzkCwPc2XIHTZiV29hAnv/0RMJIEbv4igXWfm/T6CiFEuZAwKqJ33jGvN7RkyZIR97vn8TeJJVN8fOVcrl0UIBULcerbHyEVDlK16sPMuuN/TUV1hRCibEgYFUkoFOLMmTM4HA4WLFgw7H7Pt3SwY/9ZPA4r//jBFQCcfehPiZ1+A8fs5czd+BCqgFkbhBBiOpBfvSJpaTGvst7Y2IjNNvwgxa8+9TYAX3pfE/P8boLP/yvdzz+IcriZ96fbZWYFIcQFScKoSM6cOQPAvHnzht3nhaOd/PZwO9UuG39+XROJzlO0/sg8kXX2Xf9PLgEhhLhgSRgVydmzZwGYPXv2sPv8bbpV9IX3LsLvtnP2wU2kIj34rrqVmvf+wVRUUwghypKEUZGMFkYHzvbw1Ntt+JxW/sd1TXT/7sf0vfYEFo+f2X/wXZlvTghxQZMwKpLWVvOE1+HC6P+9cAyAT62ej98S5/xP7wGg4c5/wu4fvjUlhBAXAgmjImlrawNg5syZg+4LRhL8cO8pAP7k2kbaH7+PZPAsrqZrqHnPp6a0nkIIUY4kjIokGAwC4Pf7B933wz2nCMcNblxSzxJrBx1PfgOAho9/S4ZxCyEEEkZF093dDUBNTc2g+x7aexKAjWsW0P6rr6MTMarXfgzPkjVTWkchhChXEkZFMlwYvd3Wx56T3VQ5bfy3hijBF/4VlIUZt321BLUUQojyJGFUBFprenp6AKiurs6778d7TwPwkctn0/fk/WAkqVn7MZwNS6e8nkIIUa4kjIogFAphGAZutxuHw5F336OvmmF010V2gv/5A1CK+lv/uhTVFEKIsiVhVASZLrqBraLDbX0cagvhd9u55MTPwEhQddWtOGcvK0U1hRCibEkYFUE0GgXA7Xbnlf/7W+cB+OBFtXQ/sxWA2nV/OrWVE0KICiBhVASxWAwAl8uVV/7vb54D4E7nayS7TuOYvQzvxe+f8voJIUS5kzAqgkwYOZ3ObFkoluTZlg4sCpaf/CUAtTf8sUz7I4QQQ5AwKoKhwujF410kDM37GiD25lNgsVKz9mOlqqIQQpS14S+8UyJKqZVAAPAD64EtWuuW0tZqZEOF0XMtHQB8zPoyGEm8l92MrXrwVEFCCCHKMIyA3wKLtNZBpVQA2IoZSmVr6DDqBOCq9p0A1Kz9+NRXTAghKsSYw0gp1QwEtda7hrhvJbAaaAGagJah9hujRVrrYHq7E7OVVNYyo+kyYRRLGrx0vAt/qgf36d1gtVO18kOlrKIQQpS1MYWRUmod8D1gwxD3NWF2pa3PKduulGoZT/daThABbAI2F/ocU21gy+jV0z3Ekin+0PMmdKbwLn8/VndVKasohBBlbcQBDEqpJqXUVszWTucwu23C7ErLtRXYMt5KpV/3HmDnBFpYU2bg0O59p82TYNcl9wDgu/L3S1MxIYSoECO2jNItm00ASqnhWijNDA6jPcDOzA2l1EZg8QgvlRc66de9Xym1USm1M7fVVY7i8ThAdiqgV053Y9UGSzpfBKDqig+WrG5CCFEJJjSAQSnlZ4hWU3rwAUqpJq11i9Z62xifrwlo1lrfny76KbA18zwTqetkSqVSAFjS1ybad6qb5ckW7PFeHLOW4Jg1Ug4LIYSY6Gi6AAw6zpOrCXNQw1g1AXUDbgcLDSKl1DOF7D9RWuvM65IwUhw428udyTcA8Cx731RWRQghKtJEw2jwZU0nQGu9SynlT3frgTmku+znz8mEEcCh833EjRTXqYMAeJZdV6pqCSFExSi784y01jtybo6pe2+I57i+0McopfToew37epnn4FBbH0qnuCz6OiBhJIQQY1GU6YDSx44uWLlh9HZbiCbjFN5kD7bAPOz1C0tcOyGEKH8TDaPMsZy8E1NzwqlsBx1MhkwYrUgeAcC96GqZGFUIIcZgQmGUHrjQwuBjRwHGMfCgUuW2jA639bEiab5t18KVpayWEEJUjGJ00+3CnAoo18p0+QUhr5uuPcTFRiaMripltYQQomIUEkaZmbQH2szgaYIqYhqfYsmEUQwrbb2xnJaRhJEQQozFiKPp0sd+vox5vo8f2KKUWo85Y8IOyJ7gujk9fU9motSyv+xDMWXCqDvlZG7qHFU6hLVmFjb/7BLXTAghKsNo0wEFGUMLR2u9D9hXrEpVmkwY9Wg7i4zTADjnXiKDF4QQYozkSq9F0B9GDhozYTTrolJWSQghKoqEURFku+kMO43GGQAcDRJGQggxVhJGRZDbTbcw3TKSMBJCiLGTMCqCbBgZdhYaZwFwNCwtZZWEEKKiSBgVUcRIMTvVhrZYcdQvKnV1hBCiYpTdRKmVKNMyqkoGsaCx1M5H2ewlrpUQQlQOaRkVgdYaDdQkuwBw1M0vbYWEEKLCSBgVgdYabXczM9UBSBgJIUShJIyKQGtNyllFQ6odAHvtvBLXSAghKouEURGYYeRlViaMAhJGQghRCAmjIjDDyEeDYXbT2aRlJIQQBZEwKoJMGGWOGUnLSAghCiNhVARaa1J2D7W6BwBr9cwS10gIISqLhFGRaLuT2pQZRraq+hLXRgghKouEURForXHYbbiIY1idKIen1FUSQoiKImFUBFprahzmLAwJd51cx0gIIQokYVQEWmtq7Elz2x0ocW2EEKLySBgVgdYav80wt71yvEgIIQolYVQEWmtqbHEArFV1Ja6NEEJUHgmjItBaU21NAGDzSjedEEIUSsKoCLTWeCzmMSOHr6bEtRFCiMojYVQEWmu86TByev0lro0QQlQeCaMi0FrjJQqA01td4toIIUTlkTAqAq3BSwwAp1e66YQQolASRkVgaPDpCAB2j3TTCSFEoSSMiiChwZMOI4u7qsS1EUKIyiNhVAQJrfDpMAAWlxwzEkKIQkkYFUFeGEnLSAghCiZhVARJrfq76VwSRkIIUSgJoyJIaIUbczogi9Nb4toIIUTlkTAqgqRWOHQ6jOyuEtdGCCEqj4RRESRTGjsGKRRY7aWujhBCVJyyCyOl1BallFZKdSmldiqlmkpdp9GlAEgou1xYTwghxsFW6goMYbfWurJ+0VPmtYziylHiigghRGUacxgppZqBoNZ61xD3rQRWAy1AE9Ay1H7TlSVlTpKalDASQohxGVMYKaXWAd8DNgxxXxOwRWu9Pqdsu1KqRWvdMo46NWWCD1gP3Ke1Do7jeaaMBTOMEhY5XiSEGJ+WlhZaWlpYvXo1fr85rVgwGKSlpYXOzk7WrVtX4hpOrhHDKB00m4G9QOcwu20Ctg4o2wpsYYjwGoNtmfBRSnUC2zFDqWz1t4wkjIQQQ9u1axfBYBC/309LSwuBQIDdu3ezZcsWWlpa2LJlC9u2bWP79u00NzdnH3Pfffexb98+tNaAGVobNmygpaWF7du3A2ZodXZ2snPnTrZs2UJTUxPbtm0jEAjQ0tJCR0cHW7ZsKdl7H4sRwyjdstkEoJTaPMxuzQwOoz3AzswNpdRGYPEIL7Uz062X2wrSWu9Lt8rKmlWbx4wSFmeJayLE9FDpA4EywZHR0tLCvn37uOeee/LKtm41fzqbmprYunUru3blH91obm6mqamJVatWZcuamprYu3cvtbW17Nu3L7tPxqZNm9iwYQMbN27Mlq1fv55t27bllZWbCQ1gUEr5MY8R5bWatNZBpRRKqSatdYvWetsYn28l8D2t9apRdy4jVm1eclxaRkKIoQwMGTBDZf36/E6fTPfcwP2G0tTUxJEjR/LuX716NZs2bWLz5vy2w8qVK9m5c+f0DSMgAPmtmQGaMAc1jFULOa2sdKtox7hrN0WyYSTHjIQQQ1i3bh2LF5udQ7ktmYmGQ26LCfqDa/Xq1XnldXV17Nu3b0KvNdkmGkZFvXhPukXVku7WA7Nr7+5Cn0cp9Uwx6zWaTDedjKYTQgylqamJ7du3s3nzZjZv3kxTUxObNm3K67Ybj0AgMGT5UC2scld25xlV4pDwpoXz4E2wu32lrooQ08LAYy7TQXNzM83NzbS0tLBr1y62bNnCzp072blz54iP6+wcbuzY9FKUMFJK+ctp+LXW+vpCH6OUGve//pn1dQDYnJ7xPoUQYhrbtm0bt99+O36/n6amJjZu3MjGjRtZvHhxdoTdcILBsvlpnVQTnQ4oczwor62YHtiQe/+0ptNDu7FYS1sRIURZCgaDQw5iWLduXV7LZ6hutz179kxq3crFhMIo3RpqYfCxowDmbA0XRhgZmTAqu15PIUSZuO+++waVdXZ25o2GW79+Pbt37x7y8S0t+T+nmXOLBj5f5r5KU4yJUndhTgWUa2W6/IKQDSOrhJEQYrDMAIZt27axY8cOtm3bxv333z/oRNR77rmHYDCYt9/tt98OmCPnduzYkXfS65YtW7j//vsB2LFjBxs2mPMMbNiwIdsS27x5M/fddx+7du1iw4YNZRtUaqwHCpVSXcDdWusdA8r9wPYB0wHtBDZVUssoc8xoPAdOH/y/X+GavX/L0aV38MG/fqTodRNCiHKVOUF5ohNcjzYdkB/4Mub5Qn5gixc2yaEAAAUySURBVFJqPeaMCTvSFQgqpTYrpe6hf6LULZUURBPVf8xIWkZCCDEeo00HFMScm25EWut9QHmfUTWJtGGeZ4Sc9CqEEONSdhfXq0SZlpGyymg6IYQYDwmjYjAyYSTddEIIMR4SRsWQaRnJMSMhhBgXCaNikJaREEJMiIRREWhpGQkhxIRIGBWB32kOr3e7ZNZuIYQYDwmjInjPwhoAVsyuLXFNhBCiMkkYFYN00wkhxIRIGBWBTpknvcoABiGEGB8Jo2KQWbuFEGJCJIyKoH8GBgkjIYQYDwmjIshcQkKOGQkhxPhIGBVDSq5nJIQQEzHm6xlNdxO5nlEqFiKViGJxeLE4XEWvmxBClKtiXc9IwihtImEkhBAXqmKFkXTTCSGEKDkJIyGEECUnYSSEEKLkJIyEEEKUnISREEKIkpMwEkIIUXJyluYAmWGKQgghpo60jIQQQpScnPRaBEqpZwC01teXtiYil3wv5Ue+k/JUDt+LtIyEEEKUnISREEKIkpMwEkIIUXISRkIIIUpOwkgIIUTJSRgJIYQoOQkjIYQQJSdhJIQQouTkpFchhBAlJy0jIYQQJSdhJIQQouQkjIQQQpSchJEQQoiSkzASQghRchJGQgghSk7CSAghRMlJGAkhhCg5CSMhhBAlZyt1BSqdUmolsBpoAZqAFq31rtLWavpQSjUDwaE+07F89sXaR/RLfydNwOL0eqvWeseAfeS7mWJKqXXAeqAD87vZq7XeNmCf8v1etNayjHNJfwk7B5RtB5pKXbfpsADrgC5g3Xg++2LtI0veZ9MMrMy57QeOABvluynp97Ju4P8TYC9wT6V8L9JNNzGbgK0DyrYCW0pQl2lDKdWklNqK+Y++c5jdxvLZF2sf0a9Ja70vc0NrHcT8rHI/Q/lupt6mIcp2DSgv7++l1IleyQvmX4RNA8r85sda+vpNhyX9GQ/VMhr1sy/WPrLkfS57Af+A8iZAZz5H+W5K8t1sB7YMKNuC2VVX1M98sr4XaRmNk1LKzxB/uWvzL0WUUk2lqNeFYCyffbH2maz3UInSn0tTehmSfDelobXeoLXePKC4GXgUKuN7kQEM4xeA/i9hCE2YB/dE8Y3ls6dI+8h3mENrXTtE8UrMQSYtmR8j+W5KSym1Edintb4/XVT2/2ckjMbPX+oKXMDG8tkXax8xui8D96W35bspofRIx/VgtpZy7ir770XCSAgxbum/wDtz/gIXJaTNIfY7lFJ+pdRe4G6dM+CknMkxowlK96GKEhjLZ1+sfcRg6S65TVrr9UPcJ99NCaW70bYCv80tL+fvRcJo/DL9ooHcwpwvSfqzJ89YPvti7SOGtwV4/4Ay+W7Kxy7Anz4Ztuy/F+mmGyetdVAp1cLgPtQA6YO5JajWBWGsn32x9hGDpc8Du3vggWz5bqZeuoW6F1g1zOfir4TvRVpGE7MLc0qMXCvT5WJyjeWzL9Y+Ikf6ONGW3CBSSq3LGdYr383U28PgE8Qz30fmmFF5fy+lPlmrkhfMvw4GTouxE5mupJifcRfQPJ7Pvlj7yJL32TQDG+k/32gl5lQ0W+W7Ken3cg+DT0beSc6JsOX+vaj0E4lxSk8YmOmTbcIc2y9/uU1Auv/5y5ifZzPmZ7sL8z/Ajpz9Rv3si7WPyH4vXcPc3aK1Xpyzr3w3UyzdYl3M6BOlluX3ImEkhBCi5OSYkRBCiJKTMBJCCFFyEkb/f3t1LAAAAAAwyN96GHtKIgB2MgJgJyMAdjICYCcjAHYyAmAnIwB2MgJgF6VMWKKsuR/fAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAffklEQVR4nO3de5xcZZ3n8c+vu7qq752+5ULn0knIhYCRQIOACIwKAsowKrsDjisqIzKzrI7MzEsccWWccV11151VWTUIg7cBEVGDwihKEEcQSYCQhFzIlXRufUmn77eqevaPc6rT6UvS3dXddfrU9/161atOPaeqzu/kdH/z9HOeOmXOOUREJPxyMl2AiIhMDwW+iEiWUOCLiGQJBb6ISJZQ4IuIZAkFvohIllDgi4hkCQW+hI6Z3W5mG8ys18weGNRea2bOzDoG3T4zaH3MzO43szYzO2Jmd4xjm2vNbIeZJc3sg0PWfdDMEkO2e8WQutabWZeZbTezt6f1DyAyikimCxCZAoeAfwbeARSMsH6Wcy4+QvvdwDJgETAXWG9mrzrn/n0M29wE/BD44ijrn3POXTrKugeB54Br/dsjZrbMOdc4hu2KjJl6+BI6zrlHnXM/BZrH+dKbgX9yzrU457YB9wIfBDCzT5rZ82YW8R//lZltNbN8f5v3OOd+A/SMZ4Nmthw4D/isc67bOfdjYDPw3nHWLnJaCnzJRvvNrN7M/tXMqgDMrByYh9dTT9kEnO0vfxnoBe4ys2XA/wDe75wba8CvMbMmM9tpZp9J/cfhv/8e51z7KNsVmTQKfMkmTcAFeEM25wMlwA/8dcX+feug57f6z8E5lwQ+AHwMWAd8yTn30hi3+wxwDjAbr+d+E/D3g7bbOuT5A9sVmUwKfMkazrkO59wG51zcOXcUuB24ysxKgA7/aaWDXlIKtA96/T5gPVAL3DOO7e5xzu11ziWdc5uBzwE3+Ks7hmxz2HZFJosCX7JZ6lKxOc65FuAw8MZB698IbE09MLN3AhcDv8Eb4klnu+YvbwWW+P/pjLhdkcmiwJfQMbOIfzI1F8g1s3y/7U1mtsLMcsysEvgq8LRzLjWk8l28MfpyM1sJfAR4wH/PKuDbwF/indy9zsyuHbTNqL9NA/L8beb4664xszn+8krgM8DPAJxzO4GXgc/6r3k3sBr48RT+E0m2cs7ppluobnjTK92Q2914Y+d7gU683vx3gbmDXhcD7gfagKPAHYPWPQp8c9Dja/Cmf1b6j58eYZtX+Ov+l/9+ncAevCGdvEHvVeu/vhvYAbw90/+GuoXzZs7pC1BERLKBhnRERLKEAl9EJEso8EVEsoQCX0QkSwT64mlVVVWutrY202WIiMwYGzdubHLOVY+0LtCBX1tby4YNGzJdhojIjGFm+0dbpyEdEZEsocAXEckSCnwRkSwR6DF8EZFs1t/fT319PT09w792IT8/n/nz55OXlzfm95u2wDezJcCngTLn3A2ne76ISLarr6+npKSE2tpazGyg3TlHc3Mz9fX1LF68eMzvl9aQjv+Fzw1mtmVI+9X+FzrvMrM7/QL3OOduSWd7IiLZpKenh8rKypPCHsDMqKysHLHnfyrpjuE/AFw9pJBcvC+HuAZYBdxkZqvS3I6ISFYaGvanaz+VtIZ0nHPPmFntkOYLgV3OuT1+UQ8B1wOvprMtEZGZIpl0dPbFae+J09Ebp72nn/aeOF19CTp743T3J+jqS9DV67V19Z9Y7u73nvOPf3rOpI+5T8UYfg1wYNDjeuBN/hdOfB7vy5w/5Zz7wkgvNrNbgVsBFi5cOAXliYiMrjee8IK6xwvs9l4vrL02f7nXX9fT7we6v5x6Xl+cdK8839TZy9zJ2aUB03bS1jnXDNw2huetBdYC1NXV6WL9IjJm8USS9p44rd39tPX009Ydp62n33vc3T/Q424bFM7tvYOX4/TFk5NSS1E0l+L8CCX5eZTkRyiORSiKRiiM5VIYzaUwGqEgL5eiWC4F0QiFeX57LEJhNJfls0s4tL/Z++KSEYZvJvJdJlMR+AeBBYMez/fbREROKZl0dPTFae0aObDbeuLe/SiB3tmXSLuGSI55AZ0foSTmhXXJkOBOLQ9e57V7rynOj5CbM/4x9qGO5efT3Nw87MRtapZOfn7++PYt7YqGewFYZmaL8YL+RuB943kDM7sOuO7MM8+cgvJEZCo552jt7qepo4+Wrj7auocH9qg98N70hkLMoCQWoawwj9J8/1YQoawgb1BI51ESGxTUg4M7lkd+Xs6ETohOhfnz51NfX09jY+Owdal5+OOR1lccmtmDwBVAFd53dn7WOXef/+XO/4L3JdL3O+c+P5H3r6urc7p4mkjm9cYTHOvso6m9j6bOXpo7+mjq6KW5w1tu9O+b/XXx5MRzpTgWoTQ/QmnBicA+sZxHab4X4IPXpx4XRyPkTELPeiYzs43OubqR1qU7S+emUdofBx5P571FZOo452jrifuh3UdzRy9NHb00+aHd1H4ivBs7emnviY/r/UtiEapKYpQX5o0czvkn2rz1EUr9HngkV1d8mSq6tIJIiPTFkzR29NLY3ktDWw+NHb00tPWedN/Y1kNTRx99ibGfnMzNMSqLolQWx6gqjlJVHBv+2L+vKIqSn5c7hXspExXIwNcYvsgJqd54Y3sPDe1emDe29w4sN7T3DDw+3tU/5vctjkUGQjoV3tXF3n2qvao4SmVRjLKCvKwfKgmDQAa+c+4x4LG6urqPZLoWkanW3tNPfUs3B1u6qW/por6lm/qWbo76Qd7Y3kvvGKcK5uYYVcVRqktizC7Jp7o4xuzSmP84NtBeVRyjIKpeeLYJZOCLhElbTz/1x7wwP3i82w/0E8He2n36XnlRNJfZpV6AV5fGTgR5cWygfXZpjPLC6KRMB5RwUuCLpCmeSHK4tYf9zV3sa+7k9WNd7G/u5IAf8m2nOeGZn5fD/PJCamYVML+8wFsuL2BemR/wJTGKYvpVlfQF8qdIY/gSND39CQ4c62JfsxfmXqh7y/Ut3aechliQl0tNeSrMvUAffF9ZFA3MvG8Jt0AGvsbwJRO6+uLsa+pib1Mn+5o72d/c6Yd6F0faTn0Z2jmlMRZVFLGospBFlYUsrCxiYUUhC8oLqFCgS0AEMvBFplJXX5ydRzvYfriN7Ufaea2hnT2NnRxuHT3UIzlGTXkBiyqLWFThh3pFIbVVRSwoL9QJUJkRFPgSWs456lu62eYHe+p+X3PniB/fz8s1FlYUsriqmMVVXi+9trKQRRVFnDErXx8IkhlPgS+h0NUXZ9vhNrYdbmf7Ee9+x5F2OnqHnzCN5Bhnzilm5dwSVs4rZcWcEpZUF1Ezq0ChLqEWyMDXSVs5lY7eOFsPtrLlUBtbDray5WAruxs7GOm8aXVJjJVzSzhrXunA/dLqYqIRBbtkn0AGvk7aSkpbTz9bDray9WAbmw+2suVQK3ubhg/JRHKM5XOKOfuMMs6aV8LKuaWsnFdCVXEsM4WLBFAgA1+yUzLp2NXYwYZ9LWzYf4wX97ewr7lr2PPyco0Vc0t4Q00ZZ59Rxhtqylgxt0TXbxE5DQW+ZExPf4KXDxxn4/4WNuw7xsb9LcM+pBSN5HDW3BLOqSnjnBov3JfPKdGQjMgEKPBl2sQTSTYfbOXZ3c08u7uJDftahl0jZl5ZPucvKqduUTnnL6pg5bwS8nQiVWRSBDLwddI2HJJJx/Yj7Ty7u4nndjfz/N5jw2bNrJhTwoWLK6irLaeutoKaWQUZqlYk/AIZ+DppO3M1tPfw2x2N/HZnI8/ubuZYZ99J6xdXFXHx0kouWVrJRUsqdVJVZBoFMvBl5uhPJHlxfwtP72zktzsaefVw20nr55Xlc8nSKi5ZWsnFSys5Qz14kYxR4Mu4NbT38JttDfx2RyO/39VE+6Bhmvy8HC5ZWsXly6u5bHk1tZWFuo6MSEAo8GVMjrb18O9bjvCLzYd5Yd+xk+bBnzm7mCuWV3P5imouqK3Q9EiRgFLgy6gOt3bzxOYjPLHlMBv2twyEfDQ3h8uWV/EnK2dz+fJq5pcXZrZQERkTBb6cpLWrn59vPsRPXzrIC/taBtqjkRyuWF7NO1fP460rZ1OSn5fBKkVkIgIZ+JqWOb164wnWb2/kpy8d5KntDfQlvLnxsUgOf7JiNtf6IV+sb10SmdEC+RusaZlTzznHyweO88jGen6x+TDHu7zvVc0xeMuyKt69poarzp6rkBcJEf02Z5nW7n5++tJBHvzj62w/0j7QvnJuCe85r4brz61hTml+BisUkamiwM8CzjlefL2Ff3v+AL/YfIiefm/IpqIoynvW1PDe8+dz1rzSDFcpIlNNgR9i/Ykk614+xL2/23NSb/7NZ1Zy04ULuXLVHGIRTaEUyRYK/BDq6ovzwxcO8O3f7eXg8W4Aqopj/Ke6+dx4wQIWVRZluEIRyQQFfoi0dvfzwO/38cCze2nxT8IurS7itsuXcv25NbqksEiWU+CHQEdvnH/9j73c+7s9A9eTf+OCWfz1FUu58qw55OTo0gYiosCf0br7Enz3uX1887e7B3r0Fy+p5L+97UwuXlKpa9iIyEkCGfj64NWpJZKOhzcc4CtP7qSxvReAukXl3HHVci5ZWpXh6kQkqAIZ+Prg1eg27DvG3Y9tZctB7zLEq+eXcceVy7l8ebV69CJySoEMfBnuSGsPX3hiGz97+RDgXWf+U9eexXWr5ynoRWRMFPgB19Of4L7/2Ms963fR1ZcgGsnhtsuWcNsVSymM6vCJyNgpMQLKOceTrx7ln3+xjdePdQFw9dlz+fQ7z2JBhS5HLCLjp8APoP3NnXzmZ1t5ZmcjAMvnFPPZ687mzWfqhKyITJwCP0D6E0nu/d0e/u+vX6M3nqQ0P8IdVy7n/RctIpKrD02JSHoU+AHx4ust/MOjmweuefNn557BXe9aRVVxLMOViUhYKPAzLJl0fH39Lv7Pr3fiHCysKOTz7z6HtyyrznRpIhIyCvwMau3q5xMPv8xT2xswg49evoS/edtyCqK6gqWITD4Ffoa8eqiN276/kdePdTGrMI9/+fNzuWLF7EyXJSIhpsDPgPU7Grj9By/S2ZfgnJpSvvEX52uqpYhMOQX+NPv+H/bz2XVbSSQd1597Bl9872ry8zSEIyJTL5CBH8aLpznn+J9PbOdbz+wB4GNvPZNPXLlcl0UQkWkTyMndzrnHnHO3lpWVZbqUSZFMOv7hJ1v41jN7iOQYX7phNXdctUJhLyLTKpA9/DCJJ5L8/SOv8JOXDhKL5PCt/3K+Ts6KSEYo8KdQfyLJxx58iSe2HKEwmst9N1/AxUsrM12WiGQpBf4USSQdf/vwJp7YcoSS/Ajf+fCFnLewPNNliUgWU+BPAeccd/10C+s2HaIomsv3b3kTb1wwK9NliUiWC+RJ25nMOccXntjOg398nVgkh2/ffIHCXkQCQYE/iVJTL9f6s3G+8f7zNGYvIoGhIZ1J4pzjn36+jft/v5dIjvG1m9bw1pVzMl2WiMgABf4kSCQdn/7JZh564QB5ucbX33ce7zh7bqbLEhE5iQI/Tf2JJH/78CbWbTqkefYiEmgK/DTEE0k+/tBLPL75CEXRXO774AVctERj9iISTAr8CUomHXc+upnHNx+hJBbhe3/5Js7VbBwRCTDN0pkA5xyf+/mrPLKxnoK8XB748AUKexEJPAX+BHztqV088Ow+ork53PuBOs5fVJHpkkRETkuBP07rNh3iK0/uJMfga+9bw6XLqjJdkojImCjwx2Hj/hb+7kebALjrnas09VJEZhQF/hgdbu3mo9/bQF88yfsvWsiH3lyb6ZJERMZl2mbpmFkR8P+APuBp59wPpmvb6epPJLn9316iqaOPS8+s4u7rztaXl4jIjJNWD9/M7jezBjPbMqT9ajPbYWa7zOxOv/k9wCPOuY8Af5rOdqfbl3+5g437W5hbms9Xb1pDJFd/GInIzJNucj0AXD24wcxygXuAa4BVwE1mtgqYDxzwn5ZIc7vT5tevHh24GNo9f7GGiqJopksSEZmQtALfOfcMcGxI84XALufcHudcH/AQcD1Qjxf6aW93ujR39PLJH78CwCevXqnplyIyo01F8NZwoicPXtDXAI8C7zWzbwCPjfZiM7vVzDaY2YbGxsYpKG9sUl9i0tzZx8VLKrnl0sUZq0VEZDJM20lb51wn8KExPG8tsBagrq7OTXVdo1m36RBPbPGukfOlG1aTk6OTtCIys01FD/8gsGDQ4/l+24zR0tnH3eu2AnDXu1axoKIwwxWJiKRvKgL/BWCZmS02syhwI7BuPG9gZteZ2drW1tYpKO/0vvyrHbR09XPJ0kpuvGDB6V8gIjIDpDst80HgOWCFmdWb2S3OuThwO/BLYBvwsHNu63je1zn3mHPu1rKysnTKm5BX6o/z4B9fJ5JjfO56zbcXkfBIawzfOXfTKO2PA4+n896ZkEw6PvOzrTgHH37LYs6cXZLpkkREJs2MmB45XR7fcphNB44zuyTGx962LNPliIhMqkAGfibG8OOJJF/51U4APv72ZRTH9N0wIhIugQz8TIzhP7Kxnj1NndRWFvKf63SiVkTCJ5CBP9364km++pvXAPjElcvJ07VyRCSElGx4H7I61NrD8jnFXLf6jEyXIyIyJQIZ+NM5hp9MOtY+sxuAj162VJ+oFZHQCmTgT+cY/tM7G9h5tIO5pflc90b17kUkvAIZ+NPp27/bC8Atly4mGsn6fw4RCbGsTrg9jR08u7uZ/Lwc/vxCzcwRkXDL6sD/4QveVZyvW30Gpfl5Ga5GRGRqBTLwp+OkbW88wY821gNw05sWTtl2RESCIpCBPx0nbZ/a1sCxzj5Wzi1hzYJZU7YdEZGgCGTgT4efv3IYgPeeN19XxBSRrJCVgd/VF+ep7Q0AXLt6XoarERGZHlkZ+Ou3N9Ldn2DNwlnUzCrIdDkiItMikIE/1Sdtf7H5EADvfIN69yKSPQIZ+FN50vak4RwFvohkkUAG/lR6ZmcjPf1J1iycxRkazhGRLJJ1gf/0jkYA3n7WnAxXIiIyvbIq8J1zA4F/+fLqDFcjIjK9sirwdxxt50hbD9UlMc4+ozTT5YiITKusCvzBvXt92EpEsk0gA3+qpmX+YU8zAG9ZVjWp7ysiMhMEMvCnYlpmIunYuL8FgAsXV0za+4qIzBSBDPypsPNoO+09cWpmFTCvTNMxRST7ZE3gb/B793W15RmuREQkM7In8PcdA6CuVsM5IpKdsibwX3r9OADnL1QPX0SyU1YEfltPP68f6yIayWHZnOJMlyMikhFZEfjbD7cDsGJOCXm5WbHLIiLDZEX6vXrIm8+/ap4+XSsi2SuQgT/ZH7zaeqgNgFW6nIKIZLFABv5kf/Bq2xEFvohIIAN/Mjnn2N3QCcCy2TphKyLZK/SBf6Sth+7+BBVFUWYVRjNdjohIxoQ+8Pc2er37JVVFGa5ERCSzQh/4u5u8wF+swBeRLBf6wE/18BdXK/BFJLuFP/CbOgAN6YiIhD7w61u6AVhYocAXkewW6sB3znHouBf4NbN0DXwRyW6hDvy27jidfQmKormUFkQyXY6ISEaFOvAP+r37M2YV6EvLRSTrhTrwDw0KfBGRbBfIwJ+si6cdalXgi4ikBDLwJ+viaQNDOmX5k1GWiMiMFsjAnyxN7X0AzC6NZbgSEZHMC3fgd/QCUFWswBcRCXXgN3d6gV+pwBcRCXfgp4Z0qop1WWQRkdAGvnNuoIevIR0RkRAHfltPnP6EozgWIT8vN9PliIhkXGgD/8QJWw3niIhAiAO/ucMbv9cJWxERT2gD/3iXF/jlhXkZrkREJBhCG/htPXEASvMV+CIiEObA7+4HoLRAgS8iAmEO/B4/8PN1HXwREQhz4Hf7Qzrq4YuIACEO/NbUkI7G8EVEgBAH/sCQjr7aUEQECHPgq4cvInKSaQt8M1tiZveZ2SPTsb2BaZkawxcRAcYY+GZ2v5k1mNmWIe1Xm9kOM9tlZnee6j2cc3ucc7ekU+x4tPtDOiWapSMiAsBY0/AB4OvAd1MNZpYL3ANcCdQDL5jZOiAX+MKQ13/YOdeQdrXj0NWXAKAwqsAXEYExBr5z7hkzqx3SfCGwyzm3B8DMHgKud859AXjXZBY5EV193pBOYVRXyhQRgfTG8GuAA4Me1/ttIzKzSjP7JrDGzD51iufdamYbzGxDY2PjhApLJh09/UkAXRpZRMQ3beMdzrlm4LYxPG8tsBagrq7OTWRbPXFvOCcWySE3xybyFiIioZNOD/8gsGDQ4/l+W8adGL9X715EJCWdwH8BWGZmi80sCtwIrJuMoszsOjNb29raOqHXd+uErYjIMGOdlvkg8BywwszqzewW51wcuB34JbANeNg5t3UyinLOPeacu7WsrGxCr+/u9wI/Py+0nysTERm3sc7SuWmU9seBxye1okmgHr6IyHCh7AKnxvALNIYvIjIgkIGf9hh+vzcHv0BTMkVEBgQy8NMew+/z5uBrlo6IyAmBDPx0pT5lqx6+iMgJoQz8voTXw49plo6IyIBAJmK6Y/j9cS/w83IDuXsiIhkRyERMdwy/P+FdkSGqwBcRGRDKREwN6eRFQrl7IiITEspE7NOQjojIMIFMxLTH8FMnbdXDFxEZEMhETH8MP9XD16WRRURSAhn46dKQjojIcKFMxD5/lo4CX0TkhFAmYmpIJ6oxfBGRAaFMxIHAVw9fRGRAIBMx3Vk6GsMXERkukImoWToiIpMvkIGfrtRJW43hi4icEMpE7It733ilMXwRkRNCmYipi6fpWjoiIieEMhHj/hh+JEdj+CIiKaEM/ITzevi5CnwRkQGBDPx0p2X6HXxyTIEvIpISyMBPd1pmMqkevojIUIEM/HRpSEdEZLhQBn6qh68hHRGRE0IZ+Orhi4gMF8rAT7pUDz/DhYiIBEg4A1+zdEREhgll4Cc0S0dEZJhwBr7G8EVEhglk4Kf7wSvN0hERGS6QgZ/uB6/UwxcRGS6QgZ+uEz38DBciIhIg4Qx8L+/JUeKLiAwIZeAPzNLRGL6IyIBwBr7G8EVEhgll4GuWjojIcOEMfPXwRUSGCV3gO+dOnLRV3ouIDAhd4KfC3gxMQzoiIgNCF/iaoSMiMrLQBf7ApZE1niMicpLQBb56+CIiIwtk4Kdz8TTN0BERGVkgAz+di6elvvxEHXwRkZMFMvDToU/ZioiMLHyBrzF8EZERRTJdwFSoKo5SXhjNdBkiIoESusCvLomx4a4rM12GiEjghG5IR0RERqbAFxHJEgp8EZEsocAXEckSCnwRkSyhwBcRyRIKfBGRLKHAFxHJEub8a88EkZk1Avsn+PIqoGkSy8mUsOwHaF+CSvsSPOnsxyLnXPVIKwId+Okwsw3OubpM15GusOwHaF+CSvsSPFO1HxrSERHJEgp8EZEsEebAX5vpAiZJWPYDtC9BpX0JninZj9CO4YuIyMnC3MMXEZFBFPgiIlkidIFvZleb2Q4z22Vmd2a6nrEws31mttnMXjazDX5bhZk9aWav+fflfruZ2Vf9/XvFzM7LcO33m1mDmW0Z1Dbu2s3sZv/5r5nZzQHal7vN7KB/bF42s2sHrfuUvy87zOwdg9oz+jNoZgvMbL2ZvWpmW83s4377jDsup9iXmXhc8s3sj2a2yd+Xf/TbF5vZ835dPzSzqN8e8x/v8tfXnm4fT8s5F5obkAvsBpYAUWATsCrTdY2h7n1A1ZC2LwF3+st3Al/0l68FngAMuAh4PsO1XwacB2yZaO1ABbDHvy/3l8sDsi93A383wnNX+T9fMWCx/3OXG4SfQWAecJ6/XALs9OudccflFPsyE4+LAcX+ch7wvP/v/TBwo9/+TeCv/OW/Br7pL98I/PBU+ziWGsLWw78Q2OWc2+Oc6wMeAq7PcE0TdT3wHX/5O8CfDWr/rvP8AZhlZvMyUSCAc+4Z4NiQ5vHW/g7gSefcMedcC/AkcPXUV3+yUfZlNNcDDznnep1ze4FdeD9/Gf8ZdM4dds696C+3A9uAGmbgcTnFvowmyMfFOec6/Id5/s0BbwUe8duHHpfU8XoEeJuZGaPv42mFLfBrgAODHtdz6h+OoHDAr8xso5nd6rfNcc4d9pePAHP85Zmwj+OtPej7dLs/1HF/ahiEGbIv/jDAGrze5Iw+LkP2BWbgcTGzXDN7GWjA+w90N3DcORcfoa6Bmv31rUAlaexL2AJ/prrUOXcecA3wX83sssErnfd33IycPzuTa/d9A1gKnAscBv53ZssZOzMrBn4M/I1zrm3wupl2XEbYlxl5XJxzCefcucB8vF75yuncftgC/yCwYNDj+X5boDnnDvr3DcBP8H4QjqaGavz7Bv/pM2Efx1t7YPfJOXfU/yVNAvdy4k/nQO+LmeXhBeQPnHOP+s0z8riMtC8z9bikOOeOA+uBi/GG0CIj1DVQs7++DGgmjX0JW+C/ACzzz3pH8U50rMtwTadkZkVmVpJaBq4CtuDVnZoVcTPwM395HfABf2bFRUDroD/Tg2K8tf8SuMrMyv0/za/y2zJuyPmRd+MdG/D25UZ/JsViYBnwRwLwM+iP894HbHPOfWXQqhl3XEbblxl6XKrNbJa/XABciXdOYj1wg/+0occldbxuAJ7y/zIbbR9PbzrPUk/HDW/GwU68sbFPZ7qeMdS7BO+M+yZga6pmvLG63wCvAb8GKtyJM/33+Pu3GajLcP0P4v1J3Y83lnjLRGoHPox38mkX8KEA7cv3/Fpf8X/R5g16/qf9fdkBXBOUn0HgUrzhmleAl/3btTPxuJxiX2bicVkNvOTXvAX47377ErzA3gX8CIj57fn+413++iWn28fT3XRpBRGRLBG2IR0RERmFAl9EJEso8EVEsoQCX0QkSyjwRUSyhAJfRCRLKPBFRLLE/wd5V2arnNCaXAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -233,10 +293,9 @@ } ], "source": [ - "for spec in specs:\n", - " plt.plot(ells, spec[0,0], lw=2)\n", + "\n", + "plt.plot(ells, specs[2][0,0], lw=2)\n", " \n", - "plt.plot( ells, np.sum( [s[0,0] for s in specs] , axis=0), 'k-', lw=4, label='sum' )\n", "plt.legend()\n", "plt.title('150x150')\n", "plt.yscale('log')" @@ -266,7 +325,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.8.3" } }, "nbformat": 4, diff --git a/notebooks/Basic_Example.ipynb b/notebooks/Basic_Example.ipynb index e1086bb..553f2c2 100644 --- a/notebooks/Basic_Example.ipynb +++ b/notebooks/Basic_Example.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -33,7 +33,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -49,15 +49,22 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "SED arguments: (nu, nu_0, T_CMB=2.725)\n", - "Cl arguments: (ell, ell_0)\n" + "FactorizedCrossSpectrum:\n", + "- ThermalSZ (SED):\n", + " nu: null\n", + " nu_0: null\n", + "- tSZ_150_bat (Cl):\n", + " amp: 1.0\n", + " ell: null\n", + " ell_0: null\n", + "\n" ] } ], @@ -219,7 +226,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.8.3" } }, "nbformat": 4, diff --git a/notebooks/correlated_dust_sync.ipynb b/notebooks/correlated_dust_sync.ipynb index 92ed92b..a715e03 100644 --- a/notebooks/correlated_dust_sync.ipynb +++ b/notebooks/correlated_dust_sync.ipynb @@ -21,11 +21,11 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "import pysm\n", + "#import pysm\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] @@ -39,7 +39,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -62,7 +62,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -78,7 +78,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -111,11 +111,25 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 9, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CorrelatedDustSynchrotron:\n", + "- Join (SED):\n", + " - ModifiedBlackBody: {beta: null, nu: null, nu_0: null, temp: null}\n", + " - PowerLaw: {beta: null, nu: null, nu_0: null}\n", + "- CorrelatedPowerLaws (Cl): {alpha: null, amp: null, ell: null, ell_0: null, rho: null}\n", + "\n" + ] + } + ], "source": [ - "dust_sync = fgc.CorrelatedDustSynchrotron()" + "dust_sync = fgc.CorrelatedDustSynchrotron()\n", + "print(dust_sync)" ] }, { @@ -127,7 +141,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -143,7 +157,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 11, "metadata": {}, "outputs": [ { @@ -152,7 +166,7 @@ "(7, 7, 2)" ] }, - "execution_count": 8, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -247,9 +261,16 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 18, "metadata": {}, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'sed_kwargs': {'kwseq': [{'nu': None, 'nu_0': None, 'temp': None, 'beta': None}, {'nu': None, 'beta': None, 'nu_0': None}]}, 'cl_kwargs': {'ell': array([ 9., 50.]), 'alpha': array([-0.51, -1. ]), 'ell_0': 50, 'amp': array([137. , 0.6]), 'rho': 0.17}}\n" + ] + }, { "data": { "text/plain": [ @@ -257,15 +278,21 @@ "- Join (SED):\n", " - ModifiedBlackBody: {beta: null, nu: null, nu_0: null, temp: null}\n", " - PowerLaw: {beta: null, nu: null, nu_0: null}\n", - "- CorrelatedPowerLaws (Cl): {alpha: null, amp: null, ell: null, ell_0: null, rho: null}" + "- CorrelatedPowerLaws (Cl):\n", + " alpha: [-0.5099999999999998, -1.0]\n", + " amp: [137.0, 0.6]\n", + " ell: [9.0, 50.0]\n", + " ell_0: 50\n", + " rho: 0.17" ] }, - "execution_count": 11, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ + "print(dust_sync.defaults)\n", "dust_sync" ] }, @@ -278,7 +305,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ @@ -287,7 +314,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 14, "metadata": {}, "outputs": [ { @@ -305,7 +332,7 @@ " rho: 0.17" ] }, - "execution_count": 13, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -472,7 +499,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -488,7 +515,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 16, "metadata": {}, "outputs": [ { @@ -515,7 +542,7 @@ " rho: 0.17" ] }, - "execution_count": 36, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -568,7 +595,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.7.7" } }, "nbformat": 4,