From 0fa37d5a2e206289581ee8303b0eb2a544618a4d Mon Sep 17 00:00:00 2001 From: Davide Poletti Date: Mon, 19 Apr 2021 09:32:21 +0200 Subject: [PATCH 01/10] Untested changes to diff management 1. Agnostic about the structure of kwargs (which is the same of the output derivatives of the models): It can be a nested structure. 2. To address it, _apply applies a function to the leaf of nested dictionaries and lists 3. Drafted derivatives for CorrelatedFactorizedCrossSpectrum, Join 4. Derivatives wrt to multipoles are compressed from (ell, ..., ell) to (ell, ..., 1). This saves a lot of memory and operations, but may lead to mistakes. We have to decide whether it is worthwhile 5. Relatively minor changes here and there which, however, might break some compatibility with external codes 6. Only the variables should be passed to diff (i.e., those that are None in the defaults) and returned in the output dictionary/list --- fgspectra/cross.py | 75 ++++++++++++++++++++---------- fgspectra/frequency.py | 68 +++++++++++++++++++-------- fgspectra/model.py | 17 +++++++ fgspectra/power.py | 102 ++++++++++++++++++++++++++++++++++------- 4 files changed, 202 insertions(+), 60 deletions(-) diff --git a/fgspectra/cross.py b/fgspectra/cross.py index 0acd332..53ae5f4 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') + 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->...ijl', diff, sed, cl) + res += np.einsum('pk...i,n...j,...nkl->...ijl', diff, sed, cl) + return res + + def diff_in_cl(diff): + return np.einsum('k...i,n...j,p...nkl->...ijl', diff, sed, cl) + + 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): diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index 1a94602..63c36dc 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,7 +69,7 @@ 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 @@ -89,11 +89,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) = ( + beta.reshape(-1, 1) + * (nu / nu_0)**(beta.reshape(-1, 1) - 1.) + * (_rj2cmb(nu) / _rj2cmb(nu_0)) + ) + return {'beta': res.reshape((beta.size,)+beta.shape+nu.shape)} class Synchrotron(PowerLaw): @@ -146,17 +160,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 {} + + return {'sed': np.eye(sed.size)} class ModifiedBlackBody(Model): @@ -338,3 +350,23 @@ 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. + """ + diffs = [s.diff(kw) for kw, s in zip(kwseq, 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 diffs diff --git a/fgspectra/model.py b/fgspectra/model.py index f3f5a82..21c34fe 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): + if isinstance(list_tuple_dict_or_value, (list, tuple)): + res = [_apply(func, v) for v in list_tuple_dict_or_value] + elif isinstance(list_tuple_dict_or_value, dict): + res = {k: _apply(func, v) for k, v in list_tuple_dict_or_value.items()} + else: + return func(list_tuple_dict_or_value) + return res diff --git a/fgspectra/power.py b/fgspectra/power.py index e2410e0..481ded5 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -83,6 +83,46 @@ 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] + def diff(self, **kwargs): + """ + Parameters + ---------- + ell: float or array + Multipole + alpha: float or array + Spectral index. + ell_0: float + Reference ell + amp: float or array + Amplitude, shape must be compatible with `alpha`. + + 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._cl[..., ell] / self._cl[..., ell_0, np.newaxis] + ) + res = res.reshape( + (alpha.size,) + alpha.shape + ell.shape) + + return {'amp': res} + class tSZ_150_bat(PowerSpectrumFromFile): """PowerSpectrum for Thermal Sunyaev-Zel'dovich (Dunkley et al. 2013).""" @@ -129,7 +169,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 +187,40 @@ 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] + + 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 * alpha.reshape(-1, 1) + * (ell / ell_0)**(alpha.reshape(-1, 1) - 1.) + ) + res['alpha'] = res_alpha.reshape( + (alpha.size,) + alpha.shape + ell.shape) + + return res + class FreeCls(Model): """ @@ -200,16 +268,16 @@ 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 {} + + # Convention: avoid using eye (wasteful for large number of multipoles) + return {'cls': np.ones((cls.size, 1))} class CorrelatedPowerLaws(PowerLaw): From aecfe03891921beb414e44aee6a4ea748afdf9ba Mon Sep 17 00:00:00 2001 From: beringueb Date: Mon, 19 Apr 2021 18:42:31 +0200 Subject: [PATCH 02/10] Fixes in frequency.py and power.py --- fgspectra/frequency.py | 10 +++++----- fgspectra/power.py | 15 +++++++++------ 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index 63c36dc..c7a9eaa 100644 --- a/fgspectra/frequency.py +++ b/fgspectra/frequency.py @@ -93,7 +93,7 @@ def diff(self, **kwargs): raise NotImplementedError( 'Derivatives with respect to nu and nu_0 are not implemented') - defaults = self.defaults() + defaults = self.defaults if defaults['beta'] is not None: return {} @@ -102,7 +102,7 @@ def diff(self, **kwargs): nu_0 = defaults['nu_0'] res = np.zeros((beta.size, beta.size, nu.size)) - np.einsum('bbf->bf', res) = ( + np.einsum('bbf->bf', res)[:] = ( beta.reshape(-1, 1) * (nu / nu_0)**(beta.reshape(-1, 1) - 1.) * (_rj2cmb(nu) / _rj2cmb(nu_0)) @@ -145,7 +145,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 @@ -164,10 +164,10 @@ def diff(self, nu=None, sed=None): raise NotImplementedError( 'Derivative with respect to nu does not make sense here') - defaults = self.defaults() + defaults = self.defaults if defaults['sed'] is not None: return {} - + sed = np.array(kwargs['sed']) return {'sed': np.eye(sed.size)} diff --git a/fgspectra/power.py b/fgspectra/power.py index 481ded5..53deab1 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -105,7 +105,7 @@ def diff(self, **kwargs): raise NotImplementedError( 'Derivatives with respect to ell and ell_0 are not implemented') - defaults = self.defaults() + defaults = self.defaults if defaults['amp'] is not None: return {} @@ -115,7 +115,7 @@ def diff(self, **kwargs): ell_0 = defaults['ell_0'] res = np.zeros((amp.size, amp.size, ell.size)) - np.einsum('aal->al', res) = ( + np.einsum('aal->al', res)[:] = ( self._cl[..., ell] / self._cl[..., ell_0, np.newaxis] ) res = res.reshape( @@ -191,7 +191,7 @@ def diff(self, **kwargs): raise NotImplementedError( 'Derivatives with respect to ell and ell_0 are not implemented') - defaults = self.defaults() + defaults = self.defaults res = {} alpha = defaults['alpha'] @@ -202,6 +202,8 @@ def diff(self, **kwargs): 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') @@ -212,7 +214,7 @@ def diff(self, **kwargs): ell_0 = defaults['ell_0'] res_alpha = np.zeros((alpha.size, alpha.size, ell.size)) - np.einsum('aal->al', res_alpha) = ( + np.einsum('aal->al', res_alpha)[:] = ( amp * alpha.reshape(-1, 1) * (ell / ell_0)**(alpha.reshape(-1, 1) - 1.) ) @@ -253,7 +255,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 @@ -272,9 +274,10 @@ def diff(self, ell=None, cls=None): raise NotImplementedError( 'Derivative with respect to ell does not make sense here') - defaults = self.defaults() + 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) return {'cls': np.ones((cls.size, 1))} From 07e7713cb0b85d80fc90ce15e0e139dfb9c74605 Mon Sep 17 00:00:00 2001 From: beringueb Date: Fri, 23 Apr 2021 16:20:53 +0200 Subject: [PATCH 03/10] Reverted to previous implementation of deriv wrt cls + minor fixes. --- fgspectra/cross.py | 2 +- fgspectra/frequency.py | 4 ++-- fgspectra/power.py | 5 ++++- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/fgspectra/cross.py b/fgspectra/cross.py index 53ae5f4..0a08052 100644 --- a/fgspectra/cross.py +++ b/fgspectra/cross.py @@ -181,7 +181,7 @@ def diff(self, sed_kwargs={}, cl_kwargs={}): def diff_in_sed(diff): res = np.einsum('...i,p...j->p...ij', sed, diff) - res += np.einsum('p...ij->p...ji') + res += np.einsum('p...ij->p...ji', res) return res[..., None] * cl def diff_in_cl(diff): diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index c7a9eaa..cd679f8 100644 --- a/fgspectra/frequency.py +++ b/fgspectra/frequency.py @@ -359,7 +359,7 @@ def diff(self, kwseq=None): joined. ``kwseq[i]`` is a dictionary containing the keyword arguments of the ``i``-th SED. """ - diffs = [s.diff(kw) for kw, s in zip(kwseq, self._seds)] + diffs = [s.diff(**kw) for kw, s in zip(kwseq, self._seds)] n_comp = len(diffs) for i in range(n_comp): def expand(sed): @@ -369,4 +369,4 @@ def expand(sed): diffs[i] = _apply(expand, diffs[i]) - return diffs + return {'kwseq':diffs} diff --git a/fgspectra/power.py b/fgspectra/power.py index 53deab1..aca96d2 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -280,7 +280,10 @@ def diff(self, **kwargs): cls = np.asarray(kwargs['cls']) # Convention: avoid using eye (wasteful for large number of multipoles) - return {'cls': np.ones((cls.size, 1))} + # Reverted to old method for now + #return {'cls': np.ones((cls.size, 1))} + return {'cls': np.eye(cls.size)} + class CorrelatedPowerLaws(PowerLaw): From a7ec79aba446a61340d1d05fc447b5e288f3b541 Mon Sep 17 00:00:00 2001 From: beringueb Date: Fri, 7 May 2021 17:10:14 +0200 Subject: [PATCH 04/10] Fixed Join() and CorrelatedFactorizedSpectrum(). Implemented PowerSpectraAndCovariance. --- fgspectra/cross.py | 6 +- fgspectra/frequency.py | 5 +- fgspectra/model.py | 8 +-- fgspectra/power.py | 41 +++++++++--- notebooks/ACT_Models.ipynb | 93 +++++++++++++++++++++++----- notebooks/Basic_Example.ipynb | 21 ++++--- notebooks/correlated_dust_sync.ipynb | 69 ++++++++++++++------- 7 files changed, 183 insertions(+), 60 deletions(-) diff --git a/fgspectra/cross.py b/fgspectra/cross.py index 0a08052..cb9ffe1 100644 --- a/fgspectra/cross.py +++ b/fgspectra/cross.py @@ -266,12 +266,12 @@ def diff(self, sed_kwargs={}, cl_kwargs={}): 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->...ijl', diff, sed, cl) - res += np.einsum('pk...i,n...j,...nkl->...ijl', diff, sed, cl) + 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->...ijl', diff, sed, cl) + 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)), diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index cd679f8..e7e9c12 100644 --- a/fgspectra/frequency.py +++ b/fgspectra/frequency.py @@ -359,7 +359,10 @@ def diff(self, kwseq=None): joined. ``kwseq[i]`` is a dictionary containing the keyword arguments of the ``i``-th SED. """ - diffs = [s.diff(**kw) for kw, s in zip(kwseq, self._seds)] + 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): diff --git a/fgspectra/model.py b/fgspectra/model.py index 21c34fe..f46194f 100644 --- a/fgspectra/model.py +++ b/fgspectra/model.py @@ -244,11 +244,11 @@ def diff_kwargs2array(self, diff_kwargs): return np.concatenate(res_list) -def _apply(func, list_tuple_dict_or_value): +def _apply(func, list_tuple_dict_or_value, **kwargs): if isinstance(list_tuple_dict_or_value, (list, tuple)): - res = [_apply(func, v) for v in list_tuple_dict_or_value] + 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) for k, v in list_tuple_dict_or_value.items()} + res = {k: _apply(func, v, **kwargs) for k, v in list_tuple_dict_or_value.items()} else: - return func(list_tuple_dict_or_value) + return func(list_tuple_dict_or_value, **kwargs) return res diff --git a/fgspectra/power.py b/fgspectra/power.py index aca96d2..cf9fc25 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,7 +81,8 @@ 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): """ @@ -92,7 +93,7 @@ def diff(self, **kwargs): alpha: float or array Spectral index. ell_0: float - Reference ell + Reference ells amp: float or array Amplitude, shape must be compatible with `alpha`. @@ -357,11 +358,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): @@ -453,8 +454,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:]) @@ -470,6 +474,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, From 5d0d7f77adebed45f0072d2ccbeb5a536426152a Mon Sep 17 00:00:00 2001 From: beringueb Date: Thu, 24 Jun 2021 14:33:59 +0100 Subject: [PATCH 05/10] Fix powerlaw deriv, deriv of MBB and FreeFree. --- fgspectra/frequency.py | 139 ++++++++++++++++++++++++++++++++++++++--- fgspectra/power.py | 10 +-- 2 files changed, 133 insertions(+), 16 deletions(-) diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index e7e9c12..b6457f3 100644 --- a/fgspectra/frequency.py +++ b/fgspectra/frequency.py @@ -72,6 +72,7 @@ def eval(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 @@ -103,10 +104,10 @@ def diff(self, **kwargs): res = np.zeros((beta.size, beta.size, nu.size)) np.einsum('bbf->bf', res)[:] = ( - beta.reshape(-1, 1) - * (nu / nu_0)**(beta.reshape(-1, 1) - 1.) - * (_rj2cmb(nu) / _rj2cmb(nu_0)) - ) + 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)} @@ -203,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` @@ -237,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 @@ -275,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.""" diff --git a/fgspectra/power.py b/fgspectra/power.py index cf9fc25..bff9e42 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -90,12 +90,10 @@ def diff(self, **kwargs): ---------- ell: float or array Multipole - alpha: float or array - Spectral index. ell_0: float Reference ells amp: float or array - Amplitude, shape must be compatible with `alpha`. + Amplitude, Returns ------- @@ -119,8 +117,6 @@ def diff(self, **kwargs): np.einsum('aal->al', res)[:] = ( self._cl[..., ell] / self._cl[..., ell_0, np.newaxis] ) - res = res.reshape( - (alpha.size,) + alpha.shape + ell.shape) return {'amp': res} @@ -216,8 +212,8 @@ def diff(self, **kwargs): res_alpha = np.zeros((alpha.size, alpha.size, ell.size)) np.einsum('aal->al', res_alpha)[:] = ( - amp * alpha.reshape(-1, 1) - * (ell / ell_0)**(alpha.reshape(-1, 1) - 1.) + amp * np.log(ell / ell_0) + * (ell / ell_0)**(alpha.reshape(-1, 1)) ) res['alpha'] = res_alpha.reshape( (alpha.size,) + alpha.shape + ell.shape) From cd990260eaa2e330cbbff228ce7b865ae8fac942 Mon Sep 17 00:00:00 2001 From: beringueb Date: Fri, 25 Jun 2021 00:16:34 +0100 Subject: [PATCH 06/10] Fix diff PowerSpectrumFromFile --- fgspectra/power.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/fgspectra/power.py b/fgspectra/power.py index bff9e42..fa2936e 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -114,9 +114,7 @@ def diff(self, **kwargs): ell_0 = defaults['ell_0'] res = np.zeros((amp.size, amp.size, ell.size)) - np.einsum('aal->al', res)[:] = ( - self._cl[..., ell] / self._cl[..., ell_0, np.newaxis] - ) + np.einsum('aal->al', res)[:] = self.eval(ell=ell, ell_0=ell_0) return {'amp': res} From 36920cff19fac16e47f2abcff424d7b7ebffcdff Mon Sep 17 00:00:00 2001 From: beringueb Date: Fri, 25 Jun 2021 00:18:42 +0100 Subject: [PATCH 07/10] Fix diff PowerSpectrumFromFile --- fgspectra/power.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fgspectra/power.py b/fgspectra/power.py index fa2936e..d24edc2 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -114,7 +114,7 @@ def diff(self, **kwargs): 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) + np.einsum('aal->al', res)[:] = self.eval(ell=ell, ell_0=ell_0, amp=1) return {'amp': res} From 8d0fa05c77a7dbbb3fe36edbe07dc1a2ecee084a Mon Sep 17 00:00:00 2001 From: beringueb Date: Fri, 25 Jun 2021 00:36:04 +0100 Subject: [PATCH 08/10] Fixes PowerSpectrumFromFile --- fgspectra/power.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/fgspectra/power.py b/fgspectra/power.py index d24edc2..22a97e0 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -122,17 +122,17 @@ def diff(self, **kwargs): 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): From 8953e3a618689f3cec293ac9e3a02a7f662b9e7f Mon Sep 17 00:00:00 2001 From: beringueb Date: Tue, 29 Jun 2021 10:09:29 +0100 Subject: [PATCH 09/10] Ouputs D_ell for the noise --- fgspectra/cross.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fgspectra/cross.py b/fgspectra/cross.py index cb9ffe1..4262449 100644 --- a/fgspectra/cross.py +++ b/fgspectra/cross.py @@ -361,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 From ff3b5b9e6eb27c66521173dee7ffa566c37e17db Mon Sep 17 00:00:00 2001 From: beringueb Date: Mon, 5 Jul 2021 11:56:40 +0100 Subject: [PATCH 10/10] Shape of diff of PowerSpectrumFromFile --- fgspectra/power.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fgspectra/power.py b/fgspectra/power.py index 22a97e0..37558ad 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -115,8 +115,8 @@ def diff(self, **kwargs): res = np.zeros((amp.size, amp.size, ell.size)) np.einsum('aal->al', res)[:] = self.eval(ell=ell, ell_0=ell_0, amp=1) - - return {'amp': res} + res_amp = res.reshape((amp.size,) + amp.shape + ell.shape) + return {'amp': res_amp} class tSZ_150_bat(PowerSpectrumFromFile):