From 1310818dd810d56404fc8ceec48640d93403fb59 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Fri, 5 Jun 2020 17:42:18 -0400 Subject: [PATCH 01/16] start adding model sectors --- src/pint/models/timing_model.py | 312 ++++++++++++++++++++------------ 1 file changed, 198 insertions(+), 114 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index ae5d836d1..847ecbb47 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -82,6 +82,43 @@ ] +def get_component_type(component): + """A function to identify the component object's type. + + Parameters + ---------- + component: component instance + The component object need to be inspected. + + Note + ---- + Since a component can be an inheritance from other component We inspect + all the component object bases. "inspect getmro" method returns the + base classes (including 'object') in method resolution order. The + third level of inheritance class name is what we want. + Object --> component --> TypeComponent. (i.e. DelayComponent) + This class type is in the third to the last of the getmro returned + result. + + """ + # check component type + comp_base = inspect.getmro(component.__class__) + if comp_base[-2].__name__ != "Component": + raise TypeError( + "Class '%s' is not a Component type class." + % component.__class__.__name__ + ) + elif len(comp_base) < 3: + raise TypeError( + "'%s' class is not a subclass of 'Component' class." + % component.__class__.__name__ + ) + else: + comp_type = comp_base[-3].__name__ + return comp_type + + + class TimingModel(object): """Base class for timing models and components. @@ -354,14 +391,6 @@ def d_phase_d_delay_funcs(self): Dphase_Ddelay += cp.phase_derivs_wrt_delay return Dphase_Ddelay - def get_deriv_funcs(self, component_type): - """Return dictionary of derivative functions.""" - deriv_funcs = defaultdict(list) - for cp in getattr(self, component_type + "_list"): - for k, v in cp.deriv_funcs.items(): - deriv_funcs[k] += v - return dict(deriv_funcs) - def search_cmp_attr(self, name): """Search for an attribute in all components. @@ -378,41 +407,6 @@ def search_cmp_attr(self, name): except AttributeError: continue - def get_component_type(self, component): - """A function to identify the component object's type. - - Parameters - ---------- - component: component instance - The component object need to be inspected. - - Note - ---- - Since a component can be an inheritance from other component We inspect - all the component object bases. "inspect getmro" method returns the - base classes (including 'object') in method resolution order. The - third level of inheritance class name is what we want. - Object --> component --> TypeComponent. (i.e. DelayComponent) - This class type is in the third to the last of the getmro returned - result. - - """ - # check component type - comp_base = inspect.getmro(component.__class__) - if comp_base[-2].__name__ != "Component": - raise TypeError( - "Class '%s' is not a Component type class." - % component.__class__.__name__ - ) - elif len(comp_base) < 3: - raise TypeError( - "'%s' class is not a subclass of 'Component' class." - % component.__class__.__name__ - ) - else: - comp_type = comp_base[-3].__name__ - return comp_type - def map_component(self, component): """ Get the location of component. @@ -420,7 +414,7 @@ def map_component(self, component): ---------- component: str or `Component` object Component name or component object. - + Returns ------- comp: `Component` object @@ -444,7 +438,7 @@ def map_component(self, component): ) else: comp = component - comp_type = self.get_component_type(comp) + comp_type = get_component_type(comp) host_list = getattr(self, comp_type + "_list") order = host_list.index(comp) return comp, order, host_list, comp_type @@ -462,7 +456,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr If true, add a duplicate component. Default is False. """ - comp_type = self.get_component_type(component) + comp_type = get_component_type(component) if comp_type in self.component_types: comp_list = getattr(self, comp_type + "_list") cur_cps = [] @@ -504,8 +498,8 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr self.validate() def remove_component(self, component): - """ Remove one component from the timing model. - + """ Remove one component from the timing model. + Parameters ---------- component: str or `Component` object @@ -571,7 +565,7 @@ def get_components_by_category(self): def add_param_from_top(self, param, target_component, setup=False): """ Add a parameter to a timing model component. - + Parameters ---------- param: str @@ -580,7 +574,7 @@ def add_param_from_top(self, param, target_component, setup=False): Parameter host component name. If given as "" it would add parameter to the top level `TimingModel` class setup: bool, optional - Flag to run setup() function. + Flag to run setup() function. """ if target_component == "": setattr(self, param.name, param) @@ -613,7 +607,7 @@ def remove_param(self, param): self.components[target_component].remove_param(param) def get_params_mapping(self): - """Report whick component each parameter name comes from.""" + """Report which component each parameter name comes from.""" param_mapping = {} for p in self.top_level_params: param_mapping[p] = "timing_model" @@ -658,41 +652,6 @@ def param_help(self): for par, cp in self.get_params_mapping().items() ) - def delay(self, toas, cutoff_component="", include_last=True): - """Total delay for the TOAs. - - Parameters - ---------- - toas: toa.table - The toas for analysis delays. - cutoff_component: str - The delay component name that a user wants the calculation to stop - at. - include_last: bool - If the cutoff delay component is included. - - Return the total delay which will be subtracted from the given - TOA to get time of emission at the pulsar. - - """ - delay = np.zeros(toas.ntoas) * u.second - if cutoff_component == "": - idx = len(self.DelayComponent_list) - else: - delay_names = [x.__class__.__name__ for x in self.DelayComponent_list] - if cutoff_component in delay_names: - idx = delay_names.index(cutoff_component) - if include_last: - idx += 1 - else: - raise KeyError("No delay component named '%s'." % cutoff_component) - - # Do NOT cycle through delay_funcs - cycle through components until cutoff - for dc in self.DelayComponent_list[:idx]: - for df in dc.delay_funcs_component: - delay += df(toas, delay) - return delay - def phase(self, toas, abs_phase=False): """Return the model-predicted pulse phase for the given TOAs.""" # First compute the delays to "pulsar time" @@ -848,32 +807,6 @@ def jump_flags_to_params(self, toas): getattr(self, param.name).frozen = False self.components["PhaseJump"].setup() - def get_barycentric_toas(self, toas, cutoff_component=""): - """Conveniently calculate the barycentric TOAs. - - Parameters - ---------- - toas: TOAs object - The TOAs the barycentric corrections are applied on - cutoff_delay: str, optional - The cutoff delay component name. If it is not provided, it will - search for binary delay and apply all the delay before binary. - - Return - ------ - astropy.quantity. - Barycentered TOAs. - - """ - tbl = toas.table - if cutoff_component == "": - delay_list = self.DelayComponent_list - for cp in delay_list: - if cp.category == "pulsar_system": - cutoff_component = cp.__class__.__name__ - corr = self.delay(toas, cutoff_component, False) - return tbl["tdbld"] * u.day - corr - def d_phase_d_toa(self, toas, sample_step=None): """Return the derivative of phase wrt TOA. @@ -1077,7 +1010,7 @@ def designmatrix( def compare(self, othermodel, nodmx=True): """Print comparison with another model - + Parameters ---------- othermodel @@ -1087,7 +1020,7 @@ def compare(self, othermodel, nodmx=True): Returns ------- - str + str Human readable comparison, for printing """ @@ -1750,6 +1683,157 @@ def __init__(self,): self.phase_funcs_component = [] self.phase_derivs_wrt_delay = [] +class ModelSector(object): + """ A class that groups the same type of component and provide the API for + gathering the information from each component. + + Parameter + --------- + components: list of `Component` sub-class object. + Components that are belong to the same type, for example, DelayComponent. + + Note + ---- + The order of the component in the list is the order a component get computed. + """ + _apis = ('component_names', 'component_classes') + + def __init__(self, components): + # Check if the components are the same type + self.sector_name = '' + for cp in components: + cp_type = get_component_type(cp) + if self.sector_name == '': + self.sector_name = cp_type + + if cp_type != self.sector_name: + raise ValueError("Component {} is not a {} of" + " component.".format(cp.__class__.__name__, + self.sector_name)) + self.component_list = components + + @property + def component_names(self): + return [cp.__class__.__name__ for cp in self.component_list] + + @property + def component_classes(self): + return [cp.__class__ for cp in self.component_list] + + def get_quantity_funcs(self, func_list_name): + """List of all model sector functions that contribute to the final + modeled quantity. + """ + dfs = [] + for d in self.component_list: + dfs += getattr(d, func_list_name) + return dfs + + def get_deriv_funcs(self, deriv_dict_name): + """Return dictionary of derivative functions.""" + deriv_funcs = defaultdict(list) + for cp in self.component_list: + for k, v in getattr(cp, deriv_dict_name).items(): + deriv_funcs[k] += v + return dict(deriv_funcs) + + +class DelaySector(ModelSector): + """ Class for holding all delay components and their APIs + """ + _apis = ('component_names', 'component_classes', 'delay', 'delay_funcs', + 'get_barycentric_toas', 'd_delay_d_param', 'delay_deriv_funcs') + + def __init__(self, delay_components): + super(DelaySector, self).__init__(delay_components) + + @property + def delay_funcs(self): + return self.get_quantity_funcs('delay_funcs_component') + + @property + def delay_deriv_funcs(self): + """List of derivative functions for delay components.""" + return self.get_deriv_funcs("deriv_funcs") + + def delay(self, toas, cutoff_component="", include_last=True): + """Total delay for the TOAs. + + Parameters + ---------- + toas: toa.table + The toas for analysis delays. + cutoff_component: str + The delay component name that a user wants the calculation to stop + at. + include_last: bool + If the cutoff delay component is included. + + Return the total delay which will be subtracted from the given + TOA to get time of emission at the pulsar. + + """ + delay = np.zeros(toas.ntoas) * u.second + if cutoff_component == "": + idx = len(self.component_list) + else: + delay_names = self.component_names + if cutoff_component in delay_names: + idx = delay_names.index(cutoff_component) + if include_last: + idx += 1 + else: + raise KeyError("No delay component named '%s'." % cutoff_component) + + # Do NOT cycle through delay_funcs - cycle through components until cutoff + for dc in self.component_list[:idx]: + for df in dc.delay_funcs_component: + delay += df(toas, delay) + return delay + + def get_barycentric_toas(self, toas, cutoff_component=""): + """Conveniently calculate the barycentric TOAs. + + Parameters + ---------- + toas: TOAs object + The TOAs the barycentric corrections are applied on + cutoff_delay: str, optional + The cutoff delay component name. If it is not provided, it will + search for binary delay and apply all the delay before binary. + + Return + ------ + astropy.quantity. + Barycentered TOAs. + + """ + tbl = toas.table + if cutoff_component == "": + delay_list = self.component_list + for cp in delay_list: + if cp.category == "pulsar_system": + cutoff_component = cp.__class__.__name__ + corr = self.delay(toas, cutoff_component, False) + return tbl["tdbld"] * u.day - corr + + + def d_delay_d_param(self, toas, param, acc_delay=None): + """Return the derivative of delay with respect to the parameter.""" + par = getattr(self, param) + result = np.longdouble(np.zeros(toas.ntoas) * u.s / par.units) + delay_derivs = self.delay_deriv_funcs + if param not in list(delay_derivs.keys()): + raise AttributeError( + "Derivative function for '%s' is not provided" + " or not registered. " % param + ) + for df in delay_derivs[param]: + result += df(toas, param, acc_delay).to( + result.unit, equivalencies=u.dimensionless_angles() + ) + return result + class TimingModelError(ValueError): """Generic base class for timing model errors.""" From 10d88a9403896cf09f46f44278ceaae22cd6c53c Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Sat, 6 Jun 2020 17:01:06 -0400 Subject: [PATCH 02/16] add phase sector and timing model interaction --- src/pint/models/timing_model.py | 38 +++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 847ecbb47..5889858db 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -27,7 +27,7 @@ ) -__all__ = ["DEFAULT_ORDER", "TimingModel"] +__all__ = ["DEFAULT_ORDER", "TimingModel", "Component", "ModelSector"] # Parameters or lines in parfiles we don't understand but shouldn't # complain about. These are still passed to components so that they # can use them if they want to. @@ -81,7 +81,6 @@ "wave", ] - def get_component_type(component): """A function to identify the component object's type. @@ -173,8 +172,9 @@ def __init__(self, name="", components=[]): ) self.name = name self.introduces_correlated_errors = False - self.component_types = [] self.top_level_params = [] + self.model_sectors = {} + self.add_param_from_top( strParameter( name="PSR", description="Source name", aliases=["PSRJ", "PSRB"] @@ -457,8 +457,8 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr """ comp_type = get_component_type(component) - if comp_type in self.component_types: - comp_list = getattr(self, comp_type + "_list") + if comp_type in self.model_sectors.keys(): + comp_list = self.model_sectors[comp_type].component_list cur_cps = [] for cp in comp_list: cur_cps.append((order.index(cp.category), cp)) @@ -475,7 +475,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr % component.__class__.__name__ ) else: - self.component_types.append(comp_type) + self.model_sectors[comp_type] = ModelSector([component]) cur_cps = [] # link new component to TimingModel @@ -490,7 +490,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr cur_cps.append(new_cp) cur_cps.sort(key=lambda x: x[0]) new_comp_list = [c[1] for c in cur_cps] - setattr(self, comp_type + "_list", new_comp_list) + self.model_sectors[com_type].component_list = new_comp_list # Set up components self.setup() # Validate inputs @@ -1698,7 +1698,19 @@ class ModelSector(object): """ _apis = ('component_names', 'component_classes') + def __new__(cls, component, sector_map={}): + # Map a sector subclass from component type; + all_sector_map = builtin_sector_map.update(sector_map) + try: + com_type = get_component_type(component) + cls = sector_map[com_type] + + return super().__new__(cls) + def __init__(self, components): + # If only one component is given, convert it to a list + if not isinstance(components, (list, tuple)): + components = [components,] # Check if the components are the same type self.sector_name = '' for cp in components: @@ -1835,6 +1847,18 @@ def d_delay_d_param(self, toas, param, acc_delay=None): return result +class PhaseSector(ModelSector): + """ Class for holding all phase components and their APIs + """ + def __init__(self, phase_components): + super(PhaseSector, self).__init__(phase_components) + + +builtin_sector_map = {'DelayComponent': DelaySector, + 'PhaseComponent': PhaseSector, + } + + class TimingModelError(ValueError): """Generic base class for timing model errors.""" From 847fca3d8f9430841ee886a8f531d46eaf917dc7 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Sat, 6 Jun 2020 17:01:50 -0400 Subject: [PATCH 03/16] add test --- tests/test_model_sector.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 tests/test_model_sector.py diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py new file mode 100644 index 000000000..71c2439b9 --- /dev/null +++ b/tests/test_model_sector.py @@ -0,0 +1,22 @@ +""" Various tests for ModelSector and its interaction with TimingModel +""" +from __future__ import absolute_import, division, print_function + +import pytest +import os +import numpy as np +import astropy.units as u + +from pint.models import ( + TimingModel, + ModelSector, + Component +) + + +class test_sector: + def setup(self): + self.all_components = Component.component_types + + def test_init(self): + sector = ModelSector(self.all_components['AstrometryEquatorial']) From c5ff4a374c61f3463c30c679261b98eb7ec6a13a Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Sat, 6 Jun 2020 23:20:50 -0400 Subject: [PATCH 04/16] Fix bugs --- src/pint/models/__init__.py | 2 +- src/pint/models/timing_model.py | 10 +++++++--- tests/test_model_sector.py | 2 +- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index 39abd3d23..00499a297 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -18,7 +18,7 @@ from pint.models.spindown import Spindown # Import the main timing model classes -from pint.models.timing_model import TimingModel, DEFAULT_ORDER +from pint.models.timing_model import TimingModel, DEFAULT_ORDER, ModelSector from pint.models.wave import Wave from pint.models.ifunc import IFunc diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 5889858db..98f40f1e7 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -490,7 +490,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr cur_cps.append(new_cp) cur_cps.sort(key=lambda x: x[0]) new_comp_list = [c[1] for c in cur_cps] - self.model_sectors[com_type].component_list = new_comp_list + self.model_sectors[comp_type].component_list = new_comp_list # Set up components self.setup() # Validate inputs @@ -1683,6 +1683,7 @@ def __init__(self,): self.phase_funcs_component = [] self.phase_derivs_wrt_delay = [] + class ModelSector(object): """ A class that groups the same type of component and provide the API for gathering the information from each component. @@ -1698,12 +1699,15 @@ class ModelSector(object): """ _apis = ('component_names', 'component_classes') - def __new__(cls, component, sector_map={}): + def __new__(cls, components, sector_map={}): # Map a sector subclass from component type; all_sector_map = builtin_sector_map.update(sector_map) + # Assuem the first component's type is the type for all other components. + com_type = get_component_type(components[0]) try: - com_type = get_component_type(component) cls = sector_map[com_type] + except KeyError: + ValueError("Can not find the Sector class for {}".format(com_type)) return super().__new__(cls) diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py index 71c2439b9..c8f4ecf84 100644 --- a/tests/test_model_sector.py +++ b/tests/test_model_sector.py @@ -7,7 +7,7 @@ import numpy as np import astropy.units as u -from pint.models import ( +from pint.models.timing_model import ( TimingModel, ModelSector, Component From c4bd7568b33283a50cfae008675d46001d87f264 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Sun, 7 Jun 2020 13:10:10 -0400 Subject: [PATCH 05/16] Fix __new__ and component. --- src/pint/models/timing_model.py | 33 ++++++++++++++++----------------- tests/test_model_sector.py | 11 ++++++++--- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 98f40f1e7..415d77822 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -307,19 +307,15 @@ def params_ordered(self): def components(self): """All the components indexed by name.""" comps = {} - if six.PY2: - type_list = super(TimingModel, self).__getattribute__("component_types") - else: - type_list = super().__getattribute__("component_types") - for ct in type_list: - if six.PY2: - cps_list = super(TimingModel, self).__getattribute__(ct + "_list") - else: - cps_list = super().__getattribute__(ct + "_list") - for cp in cps_list: + for k, v in self.model_sectors.items(): + for cp in v.component_list: comps[cp.__class__.__name__] = cp return comps + @property + def component_types(self): + return self.model_sector.keys() + @property def delay_funcs(self): """List of all delay functions.""" @@ -1701,17 +1697,20 @@ class ModelSector(object): def __new__(cls, components, sector_map={}): # Map a sector subclass from component type; - all_sector_map = builtin_sector_map.update(sector_map) + all_sector_map = builtin_sector_map + all_sector_map.update(sector_map) + if not isinstance(components, (list, tuple)): + components = [components,] # Assuem the first component's type is the type for all other components. com_type = get_component_type(components[0]) try: - cls = sector_map[com_type] + cls = all_sector_map[com_type] except KeyError: ValueError("Can not find the Sector class for {}".format(com_type)) return super().__new__(cls) - def __init__(self, components): + def __init__(self, components, sector_map={}): # If only one component is given, convert it to a list if not isinstance(components, (list, tuple)): components = [components,] @@ -1760,8 +1759,8 @@ class DelaySector(ModelSector): _apis = ('component_names', 'component_classes', 'delay', 'delay_funcs', 'get_barycentric_toas', 'd_delay_d_param', 'delay_deriv_funcs') - def __init__(self, delay_components): - super(DelaySector, self).__init__(delay_components) + def __init__(self, delay_components, sector_map={}): + super(DelaySector, self).__init__(delay_components, sector_map=sector_map) @property def delay_funcs(self): @@ -1854,8 +1853,8 @@ def d_delay_d_param(self, toas, param, acc_delay=None): class PhaseSector(ModelSector): """ Class for holding all phase components and their APIs """ - def __init__(self, phase_components): - super(PhaseSector, self).__init__(phase_components) + def __init__(self, phase_components, sector_map={}): + super(PhaseSector, self).__init__(phase_components, sector_map=sector_map) builtin_sector_map = {'DelayComponent': DelaySector, diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py index c8f4ecf84..a5c068bc8 100644 --- a/tests/test_model_sector.py +++ b/tests/test_model_sector.py @@ -14,9 +14,14 @@ ) -class test_sector: +class TestSector: def setup(self): self.all_components = Component.component_types - def test_init(self): - sector = ModelSector(self.all_components['AstrometryEquatorial']) + def test_sector_init(self): + sector = ModelSector(self.all_components['AstrometryEquatorial']()) + + assert sector.__class__.__name__ == 'DelaySector' + assert len(sector.component_list) == 1 + + From e686179d1ceaf4cd4db4410dc0d08a0adcd6ef10 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Mon, 8 Jun 2020 00:16:03 -0400 Subject: [PATCH 06/16] add sector timing model interaction --- src/pint/models/timing_model.py | 39 ++++++++++++++++++++++++++++----- tests/test_model_sector.py | 3 ++- 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 415d77822..39098c61a 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -238,6 +238,9 @@ def __getattr__(self, name): # AttributeError it looks like it's missing entirely errmsg = "'TimingModel' object and its component has no attribute" errmsg += " '%s'." % name + if name in self.sector_methods.keys(): + sector = self.sector_methods[name] + return getattr(self.model_sectors[sector], name) try: if six.PY2: cp = super(TimingModel, self).__getattribute__("search_cmp_attr")( @@ -251,6 +254,16 @@ def __getattr__(self, name): raise AttributeError(errmsg) except: raise AttributeError(errmsg) + + @property + def sector_methods(self): + """ Get all the sector method and reorganise it by {method_name: sector_name} + """ + md = {} + for k, v in self.model_sectors.items(): + for method in v._methods: + md[method] = k + return md @property def params(self): @@ -314,7 +327,7 @@ def components(self): @property def component_types(self): - return self.model_sector.keys() + return self.model_sectors.keys() @property def delay_funcs(self): @@ -402,7 +415,7 @@ def search_cmp_attr(self, name): return cp except AttributeError: continue - + def map_component(self, component): """ Get the location of component. @@ -438,6 +451,19 @@ def map_component(self, component): host_list = getattr(self, comp_type + "_list") order = host_list.index(comp) return comp, order, host_list, comp_type + + def _make_sector(self, component): + sector = ModelSector(component) + if sector.sector_name in self.model_sectors.keys(): + log.warn("Sector {} is already in the timing model. skip...") + return + common_method = set(sector._methods).intersection(self.sector_methods.keys()) + if len(common_method) != 0: + raise ValueError("Sector {}'s methods have the same method with the current " + "sector methods. But sector methods should be unique." + "Please check .sector_methods for the current sector" + " methods.".format(sector.sector_name)) + self.model_sectors[comp_type] = sector def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=True): """Add a component into TimingModel. @@ -471,7 +497,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr % component.__class__.__name__ ) else: - self.model_sectors[comp_type] = ModelSector([component]) + self._make_sector(component) cur_cps = [] # link new component to TimingModel @@ -1693,7 +1719,7 @@ class ModelSector(object): ---- The order of the component in the list is the order a component get computed. """ - _apis = ('component_names', 'component_classes') + _apis = tuple() def __new__(cls, components, sector_map={}): # Map a sector subclass from component type; @@ -1756,8 +1782,8 @@ def get_deriv_funcs(self, deriv_dict_name): class DelaySector(ModelSector): """ Class for holding all delay components and their APIs """ - _apis = ('component_names', 'component_classes', 'delay', 'delay_funcs', - 'get_barycentric_toas', 'd_delay_d_param', 'delay_deriv_funcs') + _methods = ('delay', 'delay_funcs', 'get_barycentric_toas', 'd_delay_d_param', + 'delay_deriv_funcs') def __init__(self, delay_components, sector_map={}): super(DelaySector, self).__init__(delay_components, sector_map=sector_map) @@ -1853,6 +1879,7 @@ def d_delay_d_param(self, toas, param, acc_delay=None): class PhaseSector(ModelSector): """ Class for holding all phase components and their APIs """ + _methods = tuple() def __init__(self, phase_components, sector_map={}): super(PhaseSector, self).__init__(phase_components, sector_map=sector_map) diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py index a5c068bc8..7041ad664 100644 --- a/tests/test_model_sector.py +++ b/tests/test_model_sector.py @@ -23,5 +23,6 @@ def test_sector_init(self): assert sector.__class__.__name__ == 'DelaySector' assert len(sector.component_list) == 1 - + assert hasattr(sector, 'delay') + assert hasattr(sector, 'delay_funcs') From 1642fa02f006f374f68fbf3803e4152f0f6d01f9 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Tue, 9 Jun 2020 23:25:41 -0400 Subject: [PATCH 07/16] add more test and fix bugs --- src/pint/models/__init__.py | 1 + src/pint/models/noise_model.py | 149 ++++++- src/pint/models/timing_model.py | 696 ++++++++++++++------------------ tests/test_model_sector.py | 13 +- tests/test_timing_model.py | 18 +- 5 files changed, 472 insertions(+), 405 deletions(-) diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index 00499a297..45ef03087 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -27,6 +27,7 @@ "StandardTimingModel", [AstrometryEquatorial(), Spindown(), DispersionDM(), SolarSystemShapiro()], ) + # BTTimingModel = generate_timing_model("BTTimingModel", # (Astrometry, Spindown, Dispersion, SolarSystemShapiro, BT)) # DDTimingModel = generate_timing_model("DDTimingModel", diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index dc2cb2fba..b753b7dad 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -7,7 +7,7 @@ from astropy import log from pint.models.parameter import floatParameter, maskParameter -from pint.models.timing_model import Component +from pint.models.timing_model import Component, ModelSector class NoiseComponent(Component): @@ -21,6 +21,119 @@ def validate(self,): super(NoiseComponent, self).validate() +class NoiseSector(ModelSector): + """ Class for holding all delay components and their APIs + """ + _methods = ('has_correlated_errors', 'covariance_matrix_funcs', 'scaled_sigma_funcs', + 'basis_funcs', 'covariance_matrix', 'scaled_sigma', 'noise_model_designmatrix', + 'noise_model_basis_weight', 'noise_model_dimensions') + + def __init__(self, noise_components): + self.sector_name = 'NoiseComponent' + super(NoiseSector, self).__init__(noise_components) + + @property + def has_correlated_errors(self): + """Whether or not this model has correlated errors.""" + if "NoiseComponent" in self.component_types: + for nc in self.component_list: + # recursive if necessary + if nc.introduces_correlated_errors: + return True + return False + + @property + def covariance_matrix_funcs(self,): + """List of covariance matrix functions.""" + return self.get_quantity_funcs('covariance_matrix_funcs') + + @property + def scaled_sigma_funcs(self,): + """List of scaled uncertainty functions.""" + return self.get_quantity_funcs('scaled_sigma_funcs') + + @property + def basis_funcs(self,): + """List of scaled uncertainty functions.""" + return self.get_quantity_funcs('basis_funcs') + + def covariance_matrix(self, toas): + """This a function to get the TOA covariance matrix for noise models. + If there is no noise model component provided, a diagonal matrix with + TOAs error as diagonal element will be returned. + """ + ntoa = toas.ntoas + tbl = toas.table + result = np.zeros((ntoa, ntoa)) + # When there is no noise model. + if len(self.covariance_matrix_funcs) == 0: + result += np.diag(tbl["error"].quantity.to(u.s).value ** 2) + return result + + for nf in self.covariance_matrix_funcs: + result += nf(toas) + return result + + def scaled_sigma(self, toas): + """This a function to get the scaled TOA uncertainties noise models. + If there is no noise model component provided, a vector with + TOAs error as values will be returned. + """ + ntoa = toas.ntoas + tbl = toas.table + result = np.zeros(ntoa) * u.us + # When there is no noise model. + if len(self.scaled_sigma_funcs) == 0: + result += tbl["error"].quantity + return result + + for nf in self.scaled_sigma_funcs: + result += nf(toas) + return result + + def noise_model_designmatrix(self, toas): + result = [] + if len(self.basis_funcs) == 0: + return None + + for nf in self.basis_funcs: + result.append(nf(toas)[0]) + return np.hstack([r for r in result]) + + def noise_model_basis_weight(self, toas): + result = [] + if len(self.basis_funcs) == 0: + return None + + for nf in self.basis_funcs: + result.append(nf(toas)[1]) + return np.hstack([r for r in result]) + + def noise_model_dimensions(self, toas): + """Returns a dictionary of correlated-noise components in the noise + model. Each entry contains a tuple (offset, size) where size is the + number of basis funtions for the component, and offset is their + starting location in the design matrix and weights vector.""" + result = {} + + # Correct results rely on this ordering being the + # same as what is done in the self.basis_funcs + # property. + if len(self.basis_funcs) > 0: + ntot = 0 + for nc in self.component_list: + bfs = nc.basis_funcs + if len(bfs) == 0: + continue + nbf = 0 + for bf in bfs: + nbf += len(bf(toas)[1]) + result[nc.category] = (ntot, nbf) + ntot += nbf + + return result + + class ScaleToaError(NoiseComponent): """Correct reported template fitting uncertainties. @@ -364,6 +477,40 @@ def pl_rn_cov_matrix(self, toas): return np.dot(Fmat * phi[None, :], Fmat.T) +class NoiseSector(ModelSector): + """ Class for holding all delay components and their APIs + """ + _methods = () + + def __init__(self, noise_components, sector_map={}): + self.sector_name = 'NoiseComponent' + super(NoiseSector, self).__init__(noise_components) + + @property + def has_correlated_errors(self): + """Whether or not this model has correlated errors.""" + if "NoiseComponent" in self.component_types: + for nc in self.NoiseComponent_list: + # recursive if necessary + if nc.introduces_correlated_errors: + return True + return False + + @property + def covariance_matrix_funcs(self,): + """List of covariance matrix functions.""" + return self.get_quantity_funcs('covariance_matrix_funcs') + + @property + def scaled_sigma_funcs(self,): + """List of scaled uncertainty functions.""" + return self.get_quantity_funcs('scaled_sigma_funcs') + + @property + def basis_funcs(self,): + """List of scaled uncertainty functions.""" + return self.get_quantity_funcs('basis_funcs') + def create_quantization_matrix(toas_table, dt=1, nmin=2): """Create quantization matrix mapping TOAs to observing epochs.""" isort = np.argsort(toas_table) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 39098c61a..fd5705cfe 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -25,9 +25,11 @@ strParameter, floatParameter, ) +from pint.models.noise_model import NoiseSector -__all__ = ["DEFAULT_ORDER", "TimingModel", "Component", "ModelSector"] +__all__ = ["DEFAULT_ORDER", "TimingModel", "Component", "ModelSector", + "builtin_sector_map"] # Parameters or lines in parfiles we don't understand but shouldn't # complain about. These are still passed to components so that they # can use them if they want to. @@ -119,7 +121,7 @@ def get_component_type(component): class TimingModel(object): - """Base class for timing models and components. + """High level interface class for timing models and components. Base-level object provides an interface for implementing pulsar timing models. A timing model contains different model components, for example @@ -134,6 +136,12 @@ class TimingModel(object): The model components for timing model. The order of the components in timing model will follow the order of input. + external_sector_map: dict, optional + The sector map for non-built-in sectors. It follows the convention, + {'Component type': ModelSector sub-class for this component}. Default is + '{}'. The built-in sectors are listed in the + `timing_model.builtin_sector_map`. + Notes ----- PINT models pulsar pulse time of arrival at observer from its emission process and @@ -165,7 +173,7 @@ class TimingModel(object): """ - def __init__(self, name="", components=[]): + def __init__(self, name="", components=[], external_sector_map={}): if not isinstance(name, str): raise ValueError( "First parameter should be the model name, was {!r}".format(name) @@ -192,7 +200,8 @@ def __init__(self, name="", components=[]): ) for cp in components: - self.add_component(cp, validate=False) + self.add_component(cp, validate=False, + external_sector_map=external_sector_map) def __repr__(self): return "{}(\n {}\n)".format( @@ -254,7 +263,7 @@ def __getattr__(self, name): raise AttributeError(errmsg) except: raise AttributeError(errmsg) - + @property def sector_methods(self): """ Get all the sector method and reorganise it by {method_name: sector_name} @@ -329,77 +338,6 @@ def components(self): def component_types(self): return self.model_sectors.keys() - @property - def delay_funcs(self): - """List of all delay functions.""" - dfs = [] - for d in self.DelayComponent_list: - dfs += d.delay_funcs_component - return dfs - - @property - def phase_funcs(self): - """List of all phase functions.""" - pfs = [] - for p in self.PhaseComponent_list: - pfs += p.phase_funcs_component - return pfs - - @property - def has_correlated_errors(self): - """Whether or not this model has correlated errors.""" - if "NoiseComponent" in self.component_types: - for nc in self.NoiseComponent_list: - # recursive if necessary - if nc.introduces_correlated_errors: - return True - return False - - @property - def covariance_matrix_funcs(self,): - """List of covariance matrix functions.""" - cvfs = [] - if "NoiseComponent" in self.component_types: - for nc in self.NoiseComponent_list: - cvfs += nc.covariance_matrix_funcs - return cvfs - - @property - def scaled_sigma_funcs(self,): - """List of scaled uncertainty functions.""" - ssfs = [] - if "NoiseComponent" in self.component_types: - for nc in self.NoiseComponent_list: - ssfs += nc.scaled_sigma_funcs - return ssfs - - @property - def basis_funcs(self,): - """List of scaled uncertainty functions.""" - bfs = [] - if "NoiseComponent" in self.component_types: - for nc in self.NoiseComponent_list: - bfs += nc.basis_funcs - return bfs - - @property - def phase_deriv_funcs(self): - """List of derivative functions for phase components.""" - return self.get_deriv_funcs("PhaseComponent") - - @property - def delay_deriv_funcs(self): - """List of derivative functions for delay components.""" - return self.get_deriv_funcs("DelayComponent") - - @property - def d_phase_d_delay_funcs(self): - """List of d_phase_d_delay functions.""" - Dphase_Ddelay = [] - for cp in self.PhaseComponent_list: - Dphase_Ddelay += cp.phase_derivs_wrt_delay - return Dphase_Ddelay - def search_cmp_attr(self, name): """Search for an attribute in all components. @@ -415,7 +353,7 @@ def search_cmp_attr(self, name): return cp except AttributeError: continue - + def map_component(self, component): """ Get the location of component. @@ -448,35 +386,53 @@ def map_component(self, component): else: comp = component comp_type = get_component_type(comp) - host_list = getattr(self, comp_type + "_list") + host_list = self.model_sectors[comp_type].component_list order = host_list.index(comp) return comp, order, host_list, comp_type - - def _make_sector(self, component): - sector = ModelSector(component) + + def _make_sector(self, component, external_sector_map={}): + # Map a sector subclass from component type; + all_sector_map = builtin_sector_map + all_sector_map.update(external_sector_map) + if not isinstance(component, (list, tuple)): + components = [component,] + # Assuem the first component's type is the type for all other components. + com_type = get_component_type(components[0]) + try: + sector_cls = all_sector_map[com_type] + except KeyError: + raise ValueError("Can not find the Sector class for" + " {}".format(com_type)) + sector = sector_cls(component) if sector.sector_name in self.model_sectors.keys(): log.warn("Sector {} is already in the timing model. skip...") - return + return common_method = set(sector._methods).intersection(self.sector_methods.keys()) if len(common_method) != 0: - raise ValueError("Sector {}'s methods have the same method with the current " + raise ValueError("Sector {}'s methods have the same method with the current " "sector methods. But sector methods should be unique." "Please check .sector_methods for the current sector" " methods.".format(sector.sector_name)) - self.model_sectors[comp_type] = sector + self.model_sectors[sector.sector_name] = sector + # Link sector to the Timing model class. + self.model_sectors[sector.sector_name]._parent = self - def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=True): + def add_component(self, component, order=DEFAULT_ORDER, force=False, + validate=True, external_sector_map={}): """Add a component into TimingModel. Parameters ---------- - component : Component + component : `pint.modle.Component` object. The component to be added to the timing model. - order : list, optional + order : list, optional. The component category order list. Default is the DEFAULT_ORDER. - force : bool, optional + force : bool, optional. If true, add a duplicate component. Default is False. - + validate: bool, optional. + If true, validate the component. Default is True. + external_sector_map: dict, optional. + Non built-in sector maps. Default is {}. """ comp_type = get_component_type(component) if comp_type in self.model_sectors.keys(): @@ -497,7 +453,8 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, validate=Tr % component.__class__.__name__ ) else: - self._make_sector(component) + self._make_sector(component, + external_sector_map=external_sector_map) cur_cps = [] # link new component to TimingModel @@ -564,18 +521,7 @@ def _locate_param_host(self, components, param): return result_comp def replicate(self, components=[], copy_component=False): - new_tm = TimingModel() - for ct in self.component_types: - comp_list = getattr(self, ct + "_list").values() - if not copy_component: - # if not copied, the components' _parent will point to the new - # TimingModel class. - new_tm.setup_components(comp_list) - else: - new_comp_list = [copy.deepcopy(c) for c in comp_list] - new_tm.setup_components(new_comp_list) - new_tm.top_level_params = self.top_level_params - return new_tm + raise NotImplementError() def get_components_by_category(self): """Return a dict of this model's component objects keyed by the category name""" @@ -588,15 +534,15 @@ def get_components_by_category(self): def add_param_from_top(self, param, target_component, setup=False): """ Add a parameter to a timing model component. - Parameters - ---------- - param: str - Parameter name - target_component: str - Parameter host component name. If given as "" it would add - parameter to the top level `TimingModel` class - setup: bool, optional - Flag to run setup() function. + Parameters + ---------- + param: str + Parameter name + target_component: str + Parameter host component name. If given as "" it would add + parameter to the top level `TimingModel` class + setup: bool, optional + Flag to run setup() function. """ if target_component == "": setattr(self, param.name, param) @@ -674,111 +620,6 @@ def param_help(self): for par, cp in self.get_params_mapping().items() ) - def phase(self, toas, abs_phase=False): - """Return the model-predicted pulse phase for the given TOAs.""" - # First compute the delays to "pulsar time" - delay = self.delay(toas) - phase = Phase(np.zeros(toas.ntoas), np.zeros(toas.ntoas)) - # Then compute the relevant pulse phases - for pf in self.phase_funcs: - phase += Phase(pf(toas, delay)) - - # If the absolute phase flag is on, use the TZR parameters to compute - # the absolute phase. - if abs_phase: - if "AbsPhase" not in list(self.components.keys()): - # if no absolute phase (TZRMJD), add the component to the model and calculate it - from pint.models import absolute_phase - - self.add_component(absolute_phase.AbsPhase()) - self.make_TZR_toa( - toas - ) # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this? - tz_toa = self.get_TZR_toa(toas) - tz_delay = self.delay(tz_toa) - tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table))) - for pf in self.phase_funcs: - tz_phase += Phase(pf(tz_toa, tz_delay)) - return phase - tz_phase - else: - return phase - - def covariance_matrix(self, toas): - """This a function to get the TOA covariance matrix for noise models. - If there is no noise model component provided, a diagonal matrix with - TOAs error as diagonal element will be returned. - """ - ntoa = toas.ntoas - tbl = toas.table - result = np.zeros((ntoa, ntoa)) - # When there is no noise model. - if len(self.covariance_matrix_funcs) == 0: - result += np.diag(tbl["error"].quantity.to(u.s).value ** 2) - return result - - for nf in self.covariance_matrix_funcs: - result += nf(toas) - return result - - def scaled_sigma(self, toas): - """This a function to get the scaled TOA uncertainties noise models. - If there is no noise model component provided, a vector with - TOAs error as values will be returned. - """ - ntoa = toas.ntoas - tbl = toas.table - result = np.zeros(ntoa) * u.us - # When there is no noise model. - if len(self.scaled_sigma_funcs) == 0: - result += tbl["error"].quantity - return result - - for nf in self.scaled_sigma_funcs: - result += nf(toas) - return result - - def noise_model_designmatrix(self, toas): - result = [] - if len(self.basis_funcs) == 0: - return None - - for nf in self.basis_funcs: - result.append(nf(toas)[0]) - return np.hstack([r for r in result]) - - def noise_model_basis_weight(self, toas): - result = [] - if len(self.basis_funcs) == 0: - return None - - for nf in self.basis_funcs: - result.append(nf(toas)[1]) - return np.hstack([r for r in result]) - - def noise_model_dimensions(self, toas): - """Returns a dictionary of correlated-noise components in the noise - model. Each entry contains a tuple (offset, size) where size is the - number of basis funtions for the component, and offset is their - starting location in the design matrix and weights vector.""" - result = {} - - # Correct results rely on this ordering being the - # same as what is done in the self.basis_funcs - # property. - if len(self.basis_funcs) > 0: - ntot = 0 - for nc in self.NoiseComponent_list: - bfs = nc.basis_funcs - if len(bfs) == 0: - continue - nbf = 0 - for bf in bfs: - nbf += len(bf(toas)[1]) - result[nc.category] = (ntot, nbf) - ntot += nbf - - return result - def jump_flags_to_params(self, toas): """convert jump flags in toas.table["flags"] to jump parameters in the model""" from . import jump @@ -829,156 +670,12 @@ def jump_flags_to_params(self, toas): getattr(self, param.name).frozen = False self.components["PhaseJump"].setup() - def d_phase_d_toa(self, toas, sample_step=None): - """Return the derivative of phase wrt TOA. - - Parameters - ---------- - toas : PINT TOAs class - The toas when the derivative of phase will be evaluated at. - sample_step : float optional - Finite difference steps. If not specified, it will take 1/10 of the - spin period. - - """ - copy_toas = copy.deepcopy(toas) - if sample_step is None: - pulse_period = 1.0 / (self.F0.quantity) - sample_step = pulse_period * 1000 - sample_dt = [-sample_step, 2 * sample_step] - - sample_phase = [] - for dt in sample_dt: - dt_array = [dt.value] * copy_toas.ntoas * dt._unit - deltaT = time.TimeDelta(dt_array) - copy_toas.adjust_TOAs(deltaT) - phase = self.phase(copy_toas) - sample_phase.append(phase) - # Use finite difference method. - # phase'(t) = (phase(t+h)-phase(t-h))/2+ 1/6*F2*h^2 + .. - # The error should be near 1/6*F2*h^2 - dp = sample_phase[1] - sample_phase[0] - d_phase_d_toa = dp.int / (2 * sample_step) + dp.frac / (2 * sample_step) - del copy_toas - return d_phase_d_toa.to(u.Hz) - def d_phase_d_tpulsar(self, toas): """Return the derivative of phase wrt time at the pulsar. NOT implemented yet. """ - pass - - def d_phase_d_param(self, toas, delay, param): - """Return the derivative of phase with respect to the parameter.""" - # TODO need to do correct chain rule stuff wrt delay derivs, etc - # Is it safe to assume that any param affecting delay only affects - # phase indirectly (and vice-versa)?? - par = getattr(self, param) - result = np.longdouble(np.zeros(toas.ntoas)) / par.units - param_phase_derivs = [] - phase_derivs = self.phase_deriv_funcs - delay_derivs = self.delay_deriv_funcs - if param in list(phase_derivs.keys()): - for df in phase_derivs[param]: - result += df(toas, param, delay).to( - result.unit, equivalencies=u.dimensionless_angles() - ) - else: - # Apply chain rule for the parameters in the delay. - # total_phase = Phase1(delay(param)) + Phase2(delay(param)) - # d_total_phase_d_param = d_Phase1/d_delay*d_delay/d_param + - # d_Phase2/d_delay*d_delay/d_param - # = (d_Phase1/d_delay + d_Phase2/d_delay) * - # d_delay_d_param - - d_delay_d_p = self.d_delay_d_param(toas, param) - dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second - for dpddf in self.d_phase_d_delay_funcs: - dpdd_result += dpddf(toas, delay) - result = dpdd_result * d_delay_d_p - return result.to(result.unit, equivalencies=u.dimensionless_angles()) - - def d_delay_d_param(self, toas, param, acc_delay=None): - """Return the derivative of delay with respect to the parameter.""" - par = getattr(self, param) - result = np.longdouble(np.zeros(toas.ntoas) << (u.s / par.units)) - delay_derivs = self.delay_deriv_funcs - if param not in list(delay_derivs.keys()): - raise AttributeError( - "Derivative function for '%s' is not provided" - " or not registered. " % param - ) - for df in delay_derivs[param]: - result += df(toas, param, acc_delay).to( - result.unit, equivalencies=u.dimensionless_angles() - ) - return result - - def d_phase_d_param_num(self, toas, param, step=1e-2): - """Return the derivative of phase with respect to the parameter. - - Compute the value numerically, using a symmetric finite difference. - - """ - # TODO : We need to know the range of parameter. - par = getattr(self, param) - ori_value = par.value - unit = par.units - if ori_value == 0: - h = 1.0 * step - else: - h = ori_value * step - parv = [par.value - h, par.value + h] - - phase_i = ( - np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled - ) - phase_f = ( - np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled - ) - for ii, val in enumerate(parv): - par.value = val - ph = self.phase(toas) - phase_i[:, ii] = ph.int - phase_f[:, ii] = ph.frac - res_i = -phase_i[:, 0] + phase_i[:, 1] - res_f = -phase_f[:, 0] + phase_f[:, 1] - result = (res_i + res_f) / (2.0 * h * unit) - # shift value back to the original value - par.quantity = ori_value - return result - - def d_delay_d_param_num(self, toas, param, step=1e-2): - """Return the derivative of delay with respect to the parameter. - - Compute the value numerically, using a symmetric finite difference. - - """ - # TODO : We need to know the range of parameter. - par = getattr(self, param) - ori_value = par.value - if ori_value is None: - # A parameter did not get to use in the model - log.warning("Parameter '%s' is not used by timing model." % param) - return np.zeros(toas.ntoas) * (u.second / par.units) - unit = par.units - if ori_value == 0: - h = 1.0 * step - else: - h = ori_value * step - parv = [par.value - h, par.value + h] - delay = np.zeros((toas.ntoas, 2)) - for ii, val in enumerate(parv): - par.value = val - try: - delay[:, ii] = self.delay(toas) - except: - par.value = ori_value - raise - d_delay = (-delay[:, 0] + delay[:, 1]) / 2.0 / h - par.value = ori_value - return d_delay * (u.second / unit) + raise NotImplementedError() def designmatrix( self, toas, acc_delay=None, scale_by_F0=True, incfrozen=False, incoffset=True @@ -1719,39 +1416,44 @@ class ModelSector(object): ---- The order of the component in the list is the order a component get computed. """ - _apis = tuple() - - def __new__(cls, components, sector_map={}): - # Map a sector subclass from component type; - all_sector_map = builtin_sector_map - all_sector_map.update(sector_map) - if not isinstance(components, (list, tuple)): - components = [components,] - # Assuem the first component's type is the type for all other components. - com_type = get_component_type(components[0]) - try: - cls = all_sector_map[com_type] - except KeyError: - ValueError("Can not find the Sector class for {}".format(com_type)) - - return super().__new__(cls) + _methods = tuple() - def __init__(self, components, sector_map={}): + def __init__(self, components): + if not hasattr(self, 'sector_name'): + raise TypeError("Please use ModelSector's subclass to " + "initialize for a specific model sector." + ) # If only one component is given, convert it to a list if not isinstance(components, (list, tuple)): components = [components,] # Check if the components are the same type - self.sector_name = '' for cp in components: cp_type = get_component_type(cp) - if self.sector_name == '': - self.sector_name = cp_type - if cp_type != self.sector_name: raise ValueError("Component {} is not a {} of" " component.".format(cp.__class__.__name__, self.sector_name)) self.component_list = components + self._parent = None + + def __getattr__(self, name): + try: + return super(ModelSector, self).__getattribute__(name) + except AttributeError: + try: + p = super(ModelSector, self).__getattribute__("_parent") + if p is None: + raise AttributeError( + "'%s' object has no attribute '%s'." + % (self.__class__.__name__, name) + ) + else: + return self._parent.__getattr__(name) + except: + raise AttributeError( + "'%s' object has no attribute '%s'." + % (self.__class__.__name__, name) + ) @property def component_names(self): @@ -1765,10 +1467,10 @@ def get_quantity_funcs(self, func_list_name): """List of all model sector functions that contribute to the final modeled quantity. """ - dfs = [] - for d in self.component_list: - dfs += getattr(d, func_list_name) - return dfs + fs = [] + for cp in self.component_list: + fs += getattr(cp, func_list_name) + return fs def get_deriv_funcs(self, deriv_dict_name): """Return dictionary of derivative functions.""" @@ -1782,11 +1484,17 @@ def get_deriv_funcs(self, deriv_dict_name): class DelaySector(ModelSector): """ Class for holding all delay components and their APIs """ - _methods = ('delay', 'delay_funcs', 'get_barycentric_toas', 'd_delay_d_param', - 'delay_deriv_funcs') + _methods = ('delay_components','delay', 'delay_funcs', + 'get_barycentric_toas', 'd_delay_d_param', + 'delay_deriv_funcs', 'd_delay_d_param_num', 'delay') def __init__(self, delay_components, sector_map={}): - super(DelaySector, self).__init__(delay_components, sector_map=sector_map) + self.sector_name = 'DelayComponent' + super(DelaySector, self).__init__(delay_components) + + @property + def delay_components(self): + return self.component_list @property def delay_funcs(self): @@ -1875,17 +1583,221 @@ def d_delay_d_param(self, toas, param, acc_delay=None): ) return result + def d_delay_d_param_num(self, toas, param, step=1e-2): + """Return the derivative of delay with respect to the parameter. + + Compute the value numerically, using a symmetric finite difference. + + """ + # TODO : We need to know the range of parameter. + par = getattr(self, param) + ori_value = par.value + if ori_value is None: + # A parameter did not get to use in the model + log.warning("Parameter '%s' is not used by timing model." % param) + return np.zeros(toas.ntoas) * (u.second / par.units) + unit = par.units + if ori_value == 0: + h = 1.0 * step + else: + h = ori_value * step + parv = [par.value - h, par.value + h] + delay = np.zeros((toas.ntoas, 2)) + for ii, val in enumerate(parv): + par.value = val + try: + delay[:, ii] = self.delay(toas) + except: + par.value = ori_value + raise + d_delay = (-delay[:, 0] + delay[:, 1]) / 2.0 / h + par.value = ori_value + return d_delay * (u.second / unit) + class PhaseSector(ModelSector): - """ Class for holding all phase components and their APIs + """ Class for holding all phase components and their APIs. + + Parameters + ---------- """ - _methods = tuple() + _methods = ('phase_components', 'phase_func', 'phase', 'phase_deriv_funcs', + 'd_phase_d_param', 'd_phase_d_param_num', 'd_phase_d_toa', + 'd_phase_d_delay_funcs', 'get_spin_frequency') + def __init__(self, phase_components, sector_map={}): - super(PhaseSector, self).__init__(phase_components, sector_map=sector_map) + self.sector_name = 'PhaseComponent' + super(PhaseSector, self).__init__(phase_components) + + @property + def phase_components(self): + return self.component_list + + @property + def phase_funcs(self): + """List of all phase functions.""" + return self.get_quantity_funcs('phase_funcs_component') + + @property + def phase_deriv_funcs(self): + """List of derivative functions for phase components.""" + return self.get_deriv_funcs("deriv_funcs") + + @property + def d_phase_d_delay_funcs(self): + """List of d_phase_d_delay functions.""" + Dphase_Ddelay = [] + for cp in self.component_list: + Dphase_Ddelay += cp.phase_derivs_wrt_delay + return Dphase_Ddelay + + def phase(self, toas, delay=None, abs_phase=False): + """Return the model-predicted pulse phase for the given TOAs. + + Parameters + ---------- + toas: `~pint.toa.TOAs` object. + TOAs for evaluating the phase. + delay: `numpy.ndarray`, optional. + Input time delay values. If not given, phase will calculate the + delay. Default is None, + abs_phase: bool + Flag for using the absolute phase. Default is False. + + Return + ------ + `~pint.phase.Phase` object. The spin phase that at given TOAs. + """ + # First compute the delays to "pulsar time" + delay = self.delay(toas) + phase = Phase(np.zeros(toas.ntoas), np.zeros(toas.ntoas)) + # Then compute the relevant pulse phases + for pf in self.phase_funcs: + phase += Phase(pf(toas, delay)) + + # If the absolute phase flag is on, use the TZR parameters to compute + # the absolute phase. + if abs_phase: + if "AbsPhase" not in list(self.components.keys()): + # if no absolute phase (TZRMJD), add the component to the model and calculate it + from pint.models import absolute_phase + + self.add_component(absolute_phase.AbsPhase()) + self.make_TZR_toa( + toas + ) # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this? + tz_toa = self.get_TZR_toa(toas) + tz_delay = self.delay(tz_toa) + tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table))) + for pf in self.phase_funcs: + tz_phase += Phase(pf(tz_toa, tz_delay)) + return phase - tz_phase + else: + return phase + + def d_phase_d_param(self, toas, delay, param): + """Return the derivative of phase with respect to the parameter.""" + # TODO need to do correct chain rule stuff wrt delay derivs, etc + # Is it safe to assume that any param affecting delay only affects + # phase indirectly (and vice-versa)?? + par = getattr(self, param) + result = np.longdouble(np.zeros(toas.ntoas)) / par.units + param_phase_derivs = [] + phase_derivs = self.phase_deriv_funcs + if param in list(phase_derivs.keys()): + for df in phase_derivs[param]: + result += df(toas, param, delay).to( + result.unit, equivalencies=u.dimensionless_angles() + ) + else: + # Apply chain rule for the parameters in the delay. + # total_phase = Phase1(delay(param)) + Phase2(delay(param)) + # d_total_phase_d_param = d_Phase1/d_delay*d_delay/d_param + + # d_Phase2/d_delay*d_delay/d_param + # = (d_Phase1/d_delay + d_Phase2/d_delay) * + # d_delay_d_param + d_delay_d_p = self.d_delay_d_param(toas, param) + dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second + for dpddf in self.d_phase_d_delay_funcs: + dpdd_result += dpddf(toas, delay) + result = dpdd_result * d_delay_d_p + return result.to(result.unit, equivalencies=u.dimensionless_angles()) + + def d_phase_d_param_num(self, toas, param, step=1e-2): + """Return the derivative of phase with respect to the parameter. + + Compute the value numerically, using a symmetric finite difference. + + """ + # TODO : We need to know the range of parameter. + par = getattr(self, param) + ori_value = par.value + unit = par.units + if ori_value == 0: + h = 1.0 * step + else: + h = ori_value * step + parv = [par.value - h, par.value + h] + + phase_i = ( + np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled + ) + phase_f = ( + np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled + ) + for ii, val in enumerate(parv): + par.value = val + ph = self.phase(toas) + phase_i[:, ii] = ph.int + phase_f[:, ii] = ph.frac + res_i = -phase_i[:, 0] + phase_i[:, 1] + res_f = -phase_f[:, 0] + phase_f[:, 1] + result = (res_i + res_f) / (2.0 * h * unit) + # shift value back to the original value + par.quantity = ori_value + return result + + def d_phase_d_toa(self, toas, sample_step=None): + """Return the derivative of phase wrt TOA. + + Parameters + ---------- + toas : PINT TOAs class + The toas when the derivative of phase will be evaluated at. + sample_step : float optional + Finite difference steps. If not specified, it will take 1/10 of the + spin period. + + """ + copy_toas = copy.deepcopy(toas) + if sample_step is None: + pulse_period = 1.0 / (self.F0.quantity) + sample_step = pulse_period * 1000 + sample_dt = [-sample_step, 2 * sample_step] + + sample_phase = [] + for dt in sample_dt: + dt_array = [dt.value] * copy_toas.ntoas * dt._unit + deltaT = time.TimeDelta(dt_array) + copy_toas.adjust_TOAs(deltaT) + phase = self.phase(copy_toas) + sample_phase.append(phase) + # Use finite difference method. + # phase'(t) = (phase(t+h)-phase(t-h))/2+ 1/6*F2*h^2 + .. + # The error should be near 1/6*F2*h^2 + dp = sample_phase[1] - sample_phase[0] + d_phase_d_toa = dp.int / (2 * sample_step) + dp.frac / (2 * sample_step) + del copy_toas + return d_phase_d_toa.to(u.Hz) + + def get_spin_frequency(self, toas=None, apparent_frequency=False): + pass + builtin_sector_map = {'DelayComponent': DelaySector, 'PhaseComponent': PhaseSector, + 'NoiseComponent': NoiseSector, } diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py index 7041ad664..42e99b36d 100644 --- a/tests/test_model_sector.py +++ b/tests/test_model_sector.py @@ -6,10 +6,12 @@ import os import numpy as np import astropy.units as u +from copy import deepcopy from pint.models.timing_model import ( TimingModel, ModelSector, + DelaySector, Component ) @@ -19,10 +21,15 @@ def setup(self): self.all_components = Component.component_types def test_sector_init(self): - sector = ModelSector(self.all_components['AstrometryEquatorial']()) - + sector = DelaySector(self.all_components['AstrometryEquatorial']()) + assert sector.__class__.__name__ == 'DelaySector' assert len(sector.component_list) == 1 assert hasattr(sector, 'delay') assert hasattr(sector, 'delay_funcs') - + + def test_copy(self): + sector = DelaySector(self.all_components['AstrometryEquatorial']()) + copy_sector = deepcopy(sector) + + assert id(sector) != id(copy_sector) diff --git a/tests/test_timing_model.py b/tests/test_timing_model.py index 9c120976a..cdcd209fe 100644 --- a/tests/test_timing_model.py +++ b/tests/test_timing_model.py @@ -27,17 +27,17 @@ def setup(self): def test_from_par(self): tm = get_model(self.parfile) assert len(tm.components) == 6 - assert len(tm.DelayComponent_list) == 4 - assert len(tm.PhaseComponent_list) == 2 + assert len(tm.delay_components) == 4 + assert len(tm.phase_components) == 2 # Check delay component order order = [] - for dcp in tm.DelayComponent_list: + for dcp in tm.delay_components: order.append(DEFAULT_ORDER.index(dcp.category)) assert all(np.diff(np.array(order)) > 0) # Check phase component order order = [] - for dcp in tm.PhaseComponent_list: + for dcp in tm.phase_components: order.append(DEFAULT_ORDER.index(dcp.category)) assert all(np.diff(np.array(order)) > 0) @@ -53,13 +53,13 @@ def test_component_input(self): # Test Delay order order = [] - for dcp in tm.DelayComponent_list: + for dcp in tm.delay_components: order.append(DEFAULT_ORDER.index(dcp.category)) assert all(np.diff(np.array(order)) > 0) # Test Phase order order = [] - for dcp in tm.PhaseComponent_list: + for dcp in tm.phase_components: order.append(DEFAULT_ORDER.index(dcp.category)) assert all(np.diff(np.array(order)) > 0) @@ -75,7 +75,7 @@ def test_add_component(self): assert cp._parent == tm # Test order - cp_pos = tm.DelayComponent_list.index(cp) + cp_pos = tm.delay_components.index(cp) assert cp_pos == 2 print(cp.params) @@ -134,7 +134,7 @@ def test_remove_component(self): # test remove by name tm.remove_component("BinaryELL1") assert not "BinaryELL1" in tm.components.keys() - assert not remove_cp in tm.DelayComponent_list + assert not remove_cp in tm.delay_components # test remove by component tm2 = TimingModel( @@ -144,4 +144,4 @@ def test_remove_component(self): remove_cp2 = tm2.components["BinaryELL1"] tm2.remove_component(remove_cp2) assert not "BinaryELL1" in tm2.components.keys() - assert not remove_cp2 in tm2.DelayComponent_list + assert not remove_cp2 in tm2.delay_components From c21f6a456509fc7efc153fdd22f4fbd19cd5bbb2 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Wed, 10 Jun 2020 17:36:57 -0400 Subject: [PATCH 08/16] Fix bugs --- src/pint/fitter.py | 14 +- src/pint/models/__init__.py | 2 +- src/pint/models/noise_model.py | 149 +-------- src/pint/models/timing_model.py | 517 +++++--------------------------- src/pint/utils.py | 37 +++ tests/test_model_manual.py | 3 +- tests/test_model_sector.py | 4 +- 7 files changed, 126 insertions(+), 600 deletions(-) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 05e4da758..f5453f9de 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -743,6 +743,7 @@ def fit_toas(self, maxiter=1, threshold=False, full_cov=False): accuracy where they both can be applied. """ chi2 = 0 + has_noise_model = 'NoiseComponent' in self.model.component_types for i in range(max(maxiter, 1)): fitp = self.get_fitparams() fitpv = self.get_fitparams_num() @@ -758,8 +759,12 @@ def fit_toas(self, maxiter=1, threshold=False, full_cov=False): # get any noise design matrices and weight vectors if not full_cov: - Mn = self.model.noise_model_designmatrix(self.toas) - phi = self.model.noise_model_basis_weight(self.toas) + if has_noise_model: + Mn = self.model.noise_model_designmatrix(self.toas) + phi = self.model.noise_model_basis_weight(self.toas) + else: + Mn = None + phi = None phiinv = np.zeros(M.shape[1]) if Mn is not None and phi is not None: phiinv = np.concatenate((phiinv, 1 / phi)) @@ -844,7 +849,10 @@ def fit_toas(self, maxiter=1, threshold=False, full_cov=False): # Compute the noise realizations if possible if not full_cov: - noise_dims = self.model.noise_model_dimensions(self.toas) + if has_noise_model: + noise_dims = self.model.noise_model_dimensions(self.toas) + else: + noise_dims = {} noise_resids = {} for comp in noise_dims.keys(): p0 = noise_dims[comp][0] + ntmpar diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index 45ef03087..88d1d60ba 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -18,7 +18,7 @@ from pint.models.spindown import Spindown # Import the main timing model classes -from pint.models.timing_model import TimingModel, DEFAULT_ORDER, ModelSector +from pint.models.timing_model import TimingModel, DEFAULT_ORDER from pint.models.wave import Wave from pint.models.ifunc import IFunc diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index b753b7dad..dc2cb2fba 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -7,7 +7,7 @@ from astropy import log from pint.models.parameter import floatParameter, maskParameter -from pint.models.timing_model import Component, ModelSector +from pint.models.timing_model import Component class NoiseComponent(Component): @@ -21,119 +21,6 @@ def validate(self,): super(NoiseComponent, self).validate() -class NoiseSector(ModelSector): - """ Class for holding all delay components and their APIs - """ - _methods = ('has_correlated_errors', 'covariance_matrix_funcs', 'scaled_sigma_funcs', - 'basis_funcs', 'covariance_matrix', 'scaled_sigma', 'noise_model_designmatrix', - 'noise_model_basis_weight', 'noise_model_dimensions') - - def __init__(self, noise_components): - self.sector_name = 'NoiseComponent' - super(NoiseSector, self).__init__(noise_components) - - @property - def has_correlated_errors(self): - """Whether or not this model has correlated errors.""" - if "NoiseComponent" in self.component_types: - for nc in self.component_list: - # recursive if necessary - if nc.introduces_correlated_errors: - return True - return False - - @property - def covariance_matrix_funcs(self,): - """List of covariance matrix functions.""" - return self.get_quantity_funcs('covariance_matrix_funcs') - - @property - def scaled_sigma_funcs(self,): - """List of scaled uncertainty functions.""" - return self.get_quantity_funcs('scaled_sigma_funcs') - - @property - def basis_funcs(self,): - """List of scaled uncertainty functions.""" - return self.get_quantity_funcs('basis_funcs') - - def covariance_matrix(self, toas): - """This a function to get the TOA covariance matrix for noise models. - If there is no noise model component provided, a diagonal matrix with - TOAs error as diagonal element will be returned. - """ - ntoa = toas.ntoas - tbl = toas.table - result = np.zeros((ntoa, ntoa)) - # When there is no noise model. - if len(self.covariance_matrix_funcs) == 0: - result += np.diag(tbl["error"].quantity.to(u.s).value ** 2) - return result - - for nf in self.covariance_matrix_funcs: - result += nf(toas) - return result - - def scaled_sigma(self, toas): - """This a function to get the scaled TOA uncertainties noise models. - If there is no noise model component provided, a vector with - TOAs error as values will be returned. - """ - ntoa = toas.ntoas - tbl = toas.table - result = np.zeros(ntoa) * u.us - # When there is no noise model. - if len(self.scaled_sigma_funcs) == 0: - result += tbl["error"].quantity - return result - - for nf in self.scaled_sigma_funcs: - result += nf(toas) - return result - - def noise_model_designmatrix(self, toas): - result = [] - if len(self.basis_funcs) == 0: - return None - - for nf in self.basis_funcs: - result.append(nf(toas)[0]) - return np.hstack([r for r in result]) - - def noise_model_basis_weight(self, toas): - result = [] - if len(self.basis_funcs) == 0: - return None - - for nf in self.basis_funcs: - result.append(nf(toas)[1]) - return np.hstack([r for r in result]) - - def noise_model_dimensions(self, toas): - """Returns a dictionary of correlated-noise components in the noise - model. Each entry contains a tuple (offset, size) where size is the - number of basis funtions for the component, and offset is their - starting location in the design matrix and weights vector.""" - result = {} - - # Correct results rely on this ordering being the - # same as what is done in the self.basis_funcs - # property. - if len(self.basis_funcs) > 0: - ntot = 0 - for nc in self.component_list: - bfs = nc.basis_funcs - if len(bfs) == 0: - continue - nbf = 0 - for bf in bfs: - nbf += len(bf(toas)[1]) - result[nc.category] = (ntot, nbf) - ntot += nbf - - return result - - class ScaleToaError(NoiseComponent): """Correct reported template fitting uncertainties. @@ -477,40 +364,6 @@ def pl_rn_cov_matrix(self, toas): return np.dot(Fmat * phi[None, :], Fmat.T) -class NoiseSector(ModelSector): - """ Class for holding all delay components and their APIs - """ - _methods = () - - def __init__(self, noise_components, sector_map={}): - self.sector_name = 'NoiseComponent' - super(NoiseSector, self).__init__(noise_components) - - @property - def has_correlated_errors(self): - """Whether or not this model has correlated errors.""" - if "NoiseComponent" in self.component_types: - for nc in self.NoiseComponent_list: - # recursive if necessary - if nc.introduces_correlated_errors: - return True - return False - - @property - def covariance_matrix_funcs(self,): - """List of covariance matrix functions.""" - return self.get_quantity_funcs('covariance_matrix_funcs') - - @property - def scaled_sigma_funcs(self,): - """List of scaled uncertainty functions.""" - return self.get_quantity_funcs('scaled_sigma_funcs') - - @property - def basis_funcs(self,): - """List of scaled uncertainty functions.""" - return self.get_quantity_funcs('basis_funcs') - def create_quantization_matrix(toas_table, dt=1, nmin=2): """Create quantization matrix mapping TOAs to observing epochs.""" isort = np.argsort(toas_table) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index fd5705cfe..c9b741843 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -6,7 +6,6 @@ import abc import copy -import inspect from collections import defaultdict import astropy.time as time @@ -18,18 +17,23 @@ import pint from pint.models.parameter import strParameter, maskParameter from pint.phase import Phase -from pint.utils import PrefixError, interesting_lines, lines_of, split_prefixed_name +from pint.utils import ( + PrefixError, + interesting_lines, + lines_of, + split_prefixed_name, + get_component_type, +) from pint.models.parameter import ( AngleParameter, prefixParameter, strParameter, floatParameter, ) -from pint.models.noise_model import NoiseSector +from pint.models.model_sector import builtin_sector_map -__all__ = ["DEFAULT_ORDER", "TimingModel", "Component", "ModelSector", - "builtin_sector_map"] +__all__ = ["DEFAULT_ORDER", "TimingModel", "Component"] # Parameters or lines in parfiles we don't understand but shouldn't # complain about. These are still passed to components so that they # can use them if they want to. @@ -83,42 +87,6 @@ "wave", ] -def get_component_type(component): - """A function to identify the component object's type. - - Parameters - ---------- - component: component instance - The component object need to be inspected. - - Note - ---- - Since a component can be an inheritance from other component We inspect - all the component object bases. "inspect getmro" method returns the - base classes (including 'object') in method resolution order. The - third level of inheritance class name is what we want. - Object --> component --> TypeComponent. (i.e. DelayComponent) - This class type is in the third to the last of the getmro returned - result. - - """ - # check component type - comp_base = inspect.getmro(component.__class__) - if comp_base[-2].__name__ != "Component": - raise TypeError( - "Class '%s' is not a Component type class." - % component.__class__.__name__ - ) - elif len(comp_base) < 3: - raise TypeError( - "'%s' class is not a subclass of 'Component' class." - % component.__class__.__name__ - ) - else: - comp_type = comp_base[-3].__name__ - return comp_type - - class TimingModel(object): """High level interface class for timing models and components. @@ -172,16 +140,16 @@ class TimingModel(object): rather than to any particular component. """ - + def __init__(self, name="", components=[], external_sector_map={}): if not isinstance(name, str): raise ValueError( "First parameter should be the model name, was {!r}".format(name) ) + self.model_sectors = {} self.name = name self.introduces_correlated_errors = False self.top_level_params = [] - self.model_sectors = {} self.add_param_from_top( strParameter( @@ -247,9 +215,15 @@ def __getattr__(self, name): # AttributeError it looks like it's missing entirely errmsg = "'TimingModel' object and its component has no attribute" errmsg += " '%s'." % name - if name in self.sector_methods.keys(): - sector = self.sector_methods[name] - return getattr(self.model_sectors[sector], name) + if six.PY2: + sector_methods = super(TimingModel, self).__getattribute__('sector_methods') + else: + sector_methods = super().__getattribute__('sector_methods') + + if name in sector_methods.keys(): + sector = sector_methods[name] + return getattr(sector, name) + try: if six.PY2: cp = super(TimingModel, self).__getattribute__("search_cmp_attr")( @@ -269,9 +243,13 @@ def sector_methods(self): """ Get all the sector method and reorganise it by {method_name: sector_name} """ md = {} - for k, v in self.model_sectors.items(): + if six.PY2: + sectors = super(TimingModel, self).__getattribute__('model_sectors') + else: + sectors = super().__getattribute__('model_sectors') + for k, v in sectors.items(): for method in v._methods: - md[method] = k + md[method] = v return md @property @@ -435,6 +413,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, Non built-in sector maps. Default is {}. """ comp_type = get_component_type(component) + print(self.model_sectors) if comp_type in self.model_sectors.keys(): comp_list = self.model_sectors[comp_type].component_list cur_cps = [] @@ -727,6 +706,50 @@ def designmatrix( M[:, mask] /= F0.value return M, params, units, scale_by_F0 + def covariance_matrix(self, toas): + """This a function to get the TOA covariance matrix for noise models. + If there is no noise model component provided, a diagonal matrix with + TOAs error as diagonal element will be returned. + """ + ntoa = toas.ntoas + tbl = toas.table + result = np.zeros((ntoa, ntoa)) + # When there is no noise model. + if 'NoiseComponent'not in self.component_types: + result += np.diag(tbl["error"].quantity.to(u.s).value ** 2) + return result + + for nf in self.covariance_matrix_funcs: + result += nf(toas) + return result + + @property + def has_correlated_errors(self): + """Whether or not this model has correlated errors.""" + if "NoiseComponent" in self.component_types: + for nc in self.model_sectors['NoiseComponent'].component_list: + # recursive if necessary + if nc.introduces_correlated_errors: + return True + return False + + def scaled_sigma(self, toas): + """This a function to get the scaled TOA uncertainties noise models. + If there is no noise model component provided, a vector with + TOAs error as values will be returned. + """ + ntoa = toas.ntoas + tbl = toas.table + result = np.zeros(ntoa) * u.us + # When there is no noise model. + if 'NoiseComponent' not in self.component_types: + result += tbl["error"].quantity + return result + + for nf in self.scaled_sigma_funcs: + result += nf(toas) + return result + def compare(self, othermodel, nodmx=True): """Print comparison with another model @@ -1403,404 +1426,6 @@ def __init__(self,): self.phase_derivs_wrt_delay = [] -class ModelSector(object): - """ A class that groups the same type of component and provide the API for - gathering the information from each component. - - Parameter - --------- - components: list of `Component` sub-class object. - Components that are belong to the same type, for example, DelayComponent. - - Note - ---- - The order of the component in the list is the order a component get computed. - """ - _methods = tuple() - - def __init__(self, components): - if not hasattr(self, 'sector_name'): - raise TypeError("Please use ModelSector's subclass to " - "initialize for a specific model sector." - ) - # If only one component is given, convert it to a list - if not isinstance(components, (list, tuple)): - components = [components,] - # Check if the components are the same type - for cp in components: - cp_type = get_component_type(cp) - if cp_type != self.sector_name: - raise ValueError("Component {} is not a {} of" - " component.".format(cp.__class__.__name__, - self.sector_name)) - self.component_list = components - self._parent = None - - def __getattr__(self, name): - try: - return super(ModelSector, self).__getattribute__(name) - except AttributeError: - try: - p = super(ModelSector, self).__getattribute__("_parent") - if p is None: - raise AttributeError( - "'%s' object has no attribute '%s'." - % (self.__class__.__name__, name) - ) - else: - return self._parent.__getattr__(name) - except: - raise AttributeError( - "'%s' object has no attribute '%s'." - % (self.__class__.__name__, name) - ) - - @property - def component_names(self): - return [cp.__class__.__name__ for cp in self.component_list] - - @property - def component_classes(self): - return [cp.__class__ for cp in self.component_list] - - def get_quantity_funcs(self, func_list_name): - """List of all model sector functions that contribute to the final - modeled quantity. - """ - fs = [] - for cp in self.component_list: - fs += getattr(cp, func_list_name) - return fs - - def get_deriv_funcs(self, deriv_dict_name): - """Return dictionary of derivative functions.""" - deriv_funcs = defaultdict(list) - for cp in self.component_list: - for k, v in getattr(cp, deriv_dict_name).items(): - deriv_funcs[k] += v - return dict(deriv_funcs) - - -class DelaySector(ModelSector): - """ Class for holding all delay components and their APIs - """ - _methods = ('delay_components','delay', 'delay_funcs', - 'get_barycentric_toas', 'd_delay_d_param', - 'delay_deriv_funcs', 'd_delay_d_param_num', 'delay') - - def __init__(self, delay_components, sector_map={}): - self.sector_name = 'DelayComponent' - super(DelaySector, self).__init__(delay_components) - - @property - def delay_components(self): - return self.component_list - - @property - def delay_funcs(self): - return self.get_quantity_funcs('delay_funcs_component') - - @property - def delay_deriv_funcs(self): - """List of derivative functions for delay components.""" - return self.get_deriv_funcs("deriv_funcs") - - def delay(self, toas, cutoff_component="", include_last=True): - """Total delay for the TOAs. - - Parameters - ---------- - toas: toa.table - The toas for analysis delays. - cutoff_component: str - The delay component name that a user wants the calculation to stop - at. - include_last: bool - If the cutoff delay component is included. - - Return the total delay which will be subtracted from the given - TOA to get time of emission at the pulsar. - - """ - delay = np.zeros(toas.ntoas) * u.second - if cutoff_component == "": - idx = len(self.component_list) - else: - delay_names = self.component_names - if cutoff_component in delay_names: - idx = delay_names.index(cutoff_component) - if include_last: - idx += 1 - else: - raise KeyError("No delay component named '%s'." % cutoff_component) - - # Do NOT cycle through delay_funcs - cycle through components until cutoff - for dc in self.component_list[:idx]: - for df in dc.delay_funcs_component: - delay += df(toas, delay) - return delay - - def get_barycentric_toas(self, toas, cutoff_component=""): - """Conveniently calculate the barycentric TOAs. - - Parameters - ---------- - toas: TOAs object - The TOAs the barycentric corrections are applied on - cutoff_delay: str, optional - The cutoff delay component name. If it is not provided, it will - search for binary delay and apply all the delay before binary. - - Return - ------ - astropy.quantity. - Barycentered TOAs. - - """ - tbl = toas.table - if cutoff_component == "": - delay_list = self.component_list - for cp in delay_list: - if cp.category == "pulsar_system": - cutoff_component = cp.__class__.__name__ - corr = self.delay(toas, cutoff_component, False) - return tbl["tdbld"] * u.day - corr - - - def d_delay_d_param(self, toas, param, acc_delay=None): - """Return the derivative of delay with respect to the parameter.""" - par = getattr(self, param) - result = np.longdouble(np.zeros(toas.ntoas) * u.s / par.units) - delay_derivs = self.delay_deriv_funcs - if param not in list(delay_derivs.keys()): - raise AttributeError( - "Derivative function for '%s' is not provided" - " or not registered. " % param - ) - for df in delay_derivs[param]: - result += df(toas, param, acc_delay).to( - result.unit, equivalencies=u.dimensionless_angles() - ) - return result - - def d_delay_d_param_num(self, toas, param, step=1e-2): - """Return the derivative of delay with respect to the parameter. - - Compute the value numerically, using a symmetric finite difference. - - """ - # TODO : We need to know the range of parameter. - par = getattr(self, param) - ori_value = par.value - if ori_value is None: - # A parameter did not get to use in the model - log.warning("Parameter '%s' is not used by timing model." % param) - return np.zeros(toas.ntoas) * (u.second / par.units) - unit = par.units - if ori_value == 0: - h = 1.0 * step - else: - h = ori_value * step - parv = [par.value - h, par.value + h] - delay = np.zeros((toas.ntoas, 2)) - for ii, val in enumerate(parv): - par.value = val - try: - delay[:, ii] = self.delay(toas) - except: - par.value = ori_value - raise - d_delay = (-delay[:, 0] + delay[:, 1]) / 2.0 / h - par.value = ori_value - return d_delay * (u.second / unit) - - -class PhaseSector(ModelSector): - """ Class for holding all phase components and their APIs. - - Parameters - ---------- - """ - _methods = ('phase_components', 'phase_func', 'phase', 'phase_deriv_funcs', - 'd_phase_d_param', 'd_phase_d_param_num', 'd_phase_d_toa', - 'd_phase_d_delay_funcs', 'get_spin_frequency') - - def __init__(self, phase_components, sector_map={}): - self.sector_name = 'PhaseComponent' - super(PhaseSector, self).__init__(phase_components) - - @property - def phase_components(self): - return self.component_list - - @property - def phase_funcs(self): - """List of all phase functions.""" - return self.get_quantity_funcs('phase_funcs_component') - - @property - def phase_deriv_funcs(self): - """List of derivative functions for phase components.""" - return self.get_deriv_funcs("deriv_funcs") - - @property - def d_phase_d_delay_funcs(self): - """List of d_phase_d_delay functions.""" - Dphase_Ddelay = [] - for cp in self.component_list: - Dphase_Ddelay += cp.phase_derivs_wrt_delay - return Dphase_Ddelay - - def phase(self, toas, delay=None, abs_phase=False): - """Return the model-predicted pulse phase for the given TOAs. - - Parameters - ---------- - toas: `~pint.toa.TOAs` object. - TOAs for evaluating the phase. - delay: `numpy.ndarray`, optional. - Input time delay values. If not given, phase will calculate the - delay. Default is None, - abs_phase: bool - Flag for using the absolute phase. Default is False. - - Return - ------ - `~pint.phase.Phase` object. The spin phase that at given TOAs. - """ - # First compute the delays to "pulsar time" - delay = self.delay(toas) - phase = Phase(np.zeros(toas.ntoas), np.zeros(toas.ntoas)) - # Then compute the relevant pulse phases - for pf in self.phase_funcs: - phase += Phase(pf(toas, delay)) - - # If the absolute phase flag is on, use the TZR parameters to compute - # the absolute phase. - if abs_phase: - if "AbsPhase" not in list(self.components.keys()): - # if no absolute phase (TZRMJD), add the component to the model and calculate it - from pint.models import absolute_phase - - self.add_component(absolute_phase.AbsPhase()) - self.make_TZR_toa( - toas - ) # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this? - tz_toa = self.get_TZR_toa(toas) - tz_delay = self.delay(tz_toa) - tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table))) - for pf in self.phase_funcs: - tz_phase += Phase(pf(tz_toa, tz_delay)) - return phase - tz_phase - else: - return phase - - def d_phase_d_param(self, toas, delay, param): - """Return the derivative of phase with respect to the parameter.""" - # TODO need to do correct chain rule stuff wrt delay derivs, etc - # Is it safe to assume that any param affecting delay only affects - # phase indirectly (and vice-versa)?? - par = getattr(self, param) - result = np.longdouble(np.zeros(toas.ntoas)) / par.units - param_phase_derivs = [] - phase_derivs = self.phase_deriv_funcs - if param in list(phase_derivs.keys()): - for df in phase_derivs[param]: - result += df(toas, param, delay).to( - result.unit, equivalencies=u.dimensionless_angles() - ) - else: - # Apply chain rule for the parameters in the delay. - # total_phase = Phase1(delay(param)) + Phase2(delay(param)) - # d_total_phase_d_param = d_Phase1/d_delay*d_delay/d_param + - # d_Phase2/d_delay*d_delay/d_param - # = (d_Phase1/d_delay + d_Phase2/d_delay) * - # d_delay_d_param - d_delay_d_p = self.d_delay_d_param(toas, param) - dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second - for dpddf in self.d_phase_d_delay_funcs: - dpdd_result += dpddf(toas, delay) - result = dpdd_result * d_delay_d_p - return result.to(result.unit, equivalencies=u.dimensionless_angles()) - - def d_phase_d_param_num(self, toas, param, step=1e-2): - """Return the derivative of phase with respect to the parameter. - - Compute the value numerically, using a symmetric finite difference. - - """ - # TODO : We need to know the range of parameter. - par = getattr(self, param) - ori_value = par.value - unit = par.units - if ori_value == 0: - h = 1.0 * step - else: - h = ori_value * step - parv = [par.value - h, par.value + h] - - phase_i = ( - np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled - ) - phase_f = ( - np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled - ) - for ii, val in enumerate(parv): - par.value = val - ph = self.phase(toas) - phase_i[:, ii] = ph.int - phase_f[:, ii] = ph.frac - res_i = -phase_i[:, 0] + phase_i[:, 1] - res_f = -phase_f[:, 0] + phase_f[:, 1] - result = (res_i + res_f) / (2.0 * h * unit) - # shift value back to the original value - par.quantity = ori_value - return result - - def d_phase_d_toa(self, toas, sample_step=None): - """Return the derivative of phase wrt TOA. - - Parameters - ---------- - toas : PINT TOAs class - The toas when the derivative of phase will be evaluated at. - sample_step : float optional - Finite difference steps. If not specified, it will take 1/10 of the - spin period. - - """ - copy_toas = copy.deepcopy(toas) - if sample_step is None: - pulse_period = 1.0 / (self.F0.quantity) - sample_step = pulse_period * 1000 - sample_dt = [-sample_step, 2 * sample_step] - - sample_phase = [] - for dt in sample_dt: - dt_array = [dt.value] * copy_toas.ntoas * dt._unit - deltaT = time.TimeDelta(dt_array) - copy_toas.adjust_TOAs(deltaT) - phase = self.phase(copy_toas) - sample_phase.append(phase) - # Use finite difference method. - # phase'(t) = (phase(t+h)-phase(t-h))/2+ 1/6*F2*h^2 + .. - # The error should be near 1/6*F2*h^2 - dp = sample_phase[1] - sample_phase[0] - d_phase_d_toa = dp.int / (2 * sample_step) + dp.frac / (2 * sample_step) - del copy_toas - return d_phase_d_toa.to(u.Hz) - - def get_spin_frequency(self, toas=None, apparent_frequency=False): - pass - - - -builtin_sector_map = {'DelayComponent': DelaySector, - 'PhaseComponent': PhaseSector, - 'NoiseComponent': NoiseSector, - } - - class TimingModelError(ValueError): """Generic base class for timing model errors.""" diff --git a/src/pint/utils.py b/src/pint/utils.py index 61c9a5f2b..bef060af6 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -9,6 +9,7 @@ import astropy.constants as const import numpy as np import six +import inspect import scipy.optimize.zeros as zeros from scipy.special import fdtrc @@ -1391,3 +1392,39 @@ def FTest(chi2_1, dof_1, chi2_2, dof_2): log.warning("Models have equal degrees of freedom, cannot preform F-test.") ft = False return ft + + +def get_component_type(component): + """A function to identify the component object's type. + + Parameters + ---------- + component: component instance + The component object need to be inspected. + + Note + ---- + Since a component can be an inheritance from other component We inspect + all the component object bases. "inspect getmro" method returns the + base classes (including 'object') in method resolution order. The + third level of inheritance class name is what we want. + Object --> component --> TypeComponent. (i.e. DelayComponent) + This class type is in the third to the last of the getmro returned + result. + + """ + # check component type + comp_base = inspect.getmro(component.__class__) + if comp_base[-2].__name__ != "Component": + raise TypeError( + "Class '%s' is not a Component type class." + % component.__class__.__name__ + ) + elif len(comp_base) < 3: + raise TypeError( + "'%s' class is not a subclass of 'Component' class." + % component.__class__.__name__ + ) + else: + comp_type = comp_base[-3].__name__ + return comp_type diff --git a/tests/test_model_manual.py b/tests/test_model_manual.py index d5e975919..1d1aedd47 100644 --- a/tests/test_model_manual.py +++ b/tests/test_model_manual.py @@ -11,6 +11,7 @@ from pint.models.spindown import Spindown from pint.models.model_builder import UnknownBinaryModel, get_model, get_model_new from pint.models.timing_model import MissingParameter, TimingModel, Component +from pint.utils import get_component_type from pinttestdata import datadir @@ -50,7 +51,7 @@ def test_component_categories(model): """ for k, v in model.components.items(): - assert model.get_component_type(v) != v.category + assert get_component_type(v) != v.category parfile = join(datadir, "J1744-1134.basic.par") diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py index 42e99b36d..d2c085d8e 100644 --- a/tests/test_model_sector.py +++ b/tests/test_model_sector.py @@ -10,9 +10,11 @@ from pint.models.timing_model import ( TimingModel, + Component, +) +from pint.models.model_sector import ( ModelSector, DelaySector, - Component ) From 8f37f0b644fb53ebf4ebbe5238fcf915639b985d Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Thu, 11 Jun 2020 09:43:18 -0400 Subject: [PATCH 09/16] add model sector file --- src/pint/models/model_sector.py | 476 ++++++++++++++++++++++++++++++++ 1 file changed, 476 insertions(+) create mode 100644 src/pint/models/model_sector.py diff --git a/src/pint/models/model_sector.py b/src/pint/models/model_sector.py new file mode 100644 index 000000000..33b7389c4 --- /dev/null +++ b/src/pint/models/model_sector.py @@ -0,0 +1,476 @@ +""" Defining Model Sector, classes for grouping the same type of model components and provide +the uniformed API +""" +import numpy as np +import astropy.units as u +from astropy import time +import copy +from collections import defaultdict + +from pint.utils import get_component_type +from pint.phase import Phase + + +class ModelSector(object): + """ A class that groups the same type of component and provide the API for + gathering the information from each component. + + Parameter + --------- + components: list of `Component` sub-class object. + Components that are belong to the same type, for example, DelayComponent. + + Note + ---- + The order of the component in the list is the order a component get computed. + """ + _methods = tuple() + + def __init__(self, components): + if not hasattr(self, 'sector_name'): + raise TypeError("Please use ModelSector's subclass to " + "initialize for a specific model sector." + ) + # If only one component is given, convert it to a list + if not isinstance(components, (list, tuple)): + components = [components,] + # Check if the components are the same type + for cp in components: + cp_type = get_component_type(cp) + if cp_type != self.sector_name: + raise ValueError("Component {} is not a {} of" + " component.".format(cp.__class__.__name__, + self.sector_name)) + self.component_list = components + self._parent = None + + def __getattr__(self, name): + try: + return super(ModelSector, self).__getattribute__(name) + except AttributeError: + try: + p = super(ModelSector, self).__getattribute__("_parent") + if p is None: + raise AttributeError( + "'%s' object has no attribute '%s'." + % (self.__class__.__name__, name) + ) + else: + return self._parent.__getattr__(name) + except: + raise AttributeError( + "'%s' object has no attribute '%s'." + % (self.__class__.__name__, name) + ) + + @property + def component_names(self): + return [cp.__class__.__name__ for cp in self.component_list] + + @property + def component_classes(self): + return [cp.__class__ for cp in self.component_list] + + def get_quantity_funcs(self, func_list_name): + """List of all model sector functions that contribute to the final + modeled quantity. + """ + fs = [] + for cp in self.component_list: + fs += getattr(cp, func_list_name) + return fs + + def get_deriv_funcs(self, deriv_dict_name): + """Return dictionary of derivative functions.""" + deriv_funcs = defaultdict(list) + for cp in self.component_list: + for k, v in getattr(cp, deriv_dict_name).items(): + deriv_funcs[k] += v + return dict(deriv_funcs) + + +class DelaySector(ModelSector): + """ Class for holding all delay components and their APIs + """ + _methods = ('delay_components','delay', 'delay_funcs', + 'get_barycentric_toas', 'd_delay_d_param', + 'delay_deriv_funcs', 'd_delay_d_param_num', 'delay') + + def __init__(self, delay_components, sector_map={}): + self.sector_name = 'DelayComponent' + super(DelaySector, self).__init__(delay_components) + + @property + def delay_components(self): + return self.component_list + + @property + def delay_funcs(self): + return self.get_quantity_funcs('delay_funcs_component') + + @property + def delay_deriv_funcs(self): + """List of derivative functions for delay components.""" + return self.get_deriv_funcs("deriv_funcs") + + def delay(self, toas, cutoff_component="", include_last=True): + """Total delay for the TOAs. + + Parameters + ---------- + toas: toa.table + The toas for analysis delays. + cutoff_component: str + The delay component name that a user wants the calculation to stop + at. + include_last: bool + If the cutoff delay component is included. + + Return the total delay which will be subtracted from the given + TOA to get time of emission at the pulsar. + + """ + delay = np.zeros(toas.ntoas) * u.second + if cutoff_component == "": + idx = len(self.component_list) + else: + delay_names = self.component_names + if cutoff_component in delay_names: + idx = delay_names.index(cutoff_component) + if include_last: + idx += 1 + else: + raise KeyError("No delay component named '%s'." % cutoff_component) + + # Do NOT cycle through delay_funcs - cycle through components until cutoff + for dc in self.component_list[:idx]: + for df in dc.delay_funcs_component: + delay += df(toas, delay) + return delay + + def get_barycentric_toas(self, toas, cutoff_component=""): + """Conveniently calculate the barycentric TOAs. + + Parameters + ---------- + toas: TOAs object + The TOAs the barycentric corrections are applied on + cutoff_delay: str, optional + The cutoff delay component name. If it is not provided, it will + search for binary delay and apply all the delay before binary. + + Return + ------ + astropy.quantity. + Barycentered TOAs. + + """ + tbl = toas.table + if cutoff_component == "": + delay_list = self.component_list + for cp in delay_list: + if cp.category == "pulsar_system": + cutoff_component = cp.__class__.__name__ + corr = self.delay(toas, cutoff_component, False) + return tbl["tdbld"] * u.day - corr + + + def d_delay_d_param(self, toas, param, acc_delay=None): + """Return the derivative of delay with respect to the parameter.""" + par = getattr(self, param) + result = np.longdouble(np.zeros(toas.ntoas) * u.s / par.units) + delay_derivs = self.delay_deriv_funcs + if param not in list(delay_derivs.keys()): + raise AttributeError( + "Derivative function for '%s' is not provided" + " or not registered. " % param + ) + for df in delay_derivs[param]: + result += df(toas, param, acc_delay).to( + result.unit, equivalencies=u.dimensionless_angles() + ) + return result + + def d_delay_d_param_num(self, toas, param, step=1e-2): + """Return the derivative of delay with respect to the parameter. + + Compute the value numerically, using a symmetric finite difference. + + """ + # TODO : We need to know the range of parameter. + par = getattr(self, param) + ori_value = par.value + if ori_value is None: + # A parameter did not get to use in the model + log.warning("Parameter '%s' is not used by timing model." % param) + return np.zeros(toas.ntoas) * (u.second / par.units) + unit = par.units + if ori_value == 0: + h = 1.0 * step + else: + h = ori_value * step + parv = [par.value - h, par.value + h] + delay = np.zeros((toas.ntoas, 2)) + for ii, val in enumerate(parv): + par.value = val + try: + delay[:, ii] = self.delay(toas) + except: + par.value = ori_value + raise + d_delay = (-delay[:, 0] + delay[:, 1]) / 2.0 / h + par.value = ori_value + return d_delay * (u.second / unit) + + +class PhaseSector(ModelSector): + """ Class for holding all phase components and their APIs. + + Parameters + ---------- + """ + _methods = ('phase_components', 'phase_func', 'phase', 'phase_deriv_funcs', + 'd_phase_d_param', 'd_phase_d_param_num', 'd_phase_d_toa', + 'd_phase_d_delay_funcs', 'get_spin_frequency') + + def __init__(self, phase_components, sector_map={}): + self.sector_name = 'PhaseComponent' + super(PhaseSector, self).__init__(phase_components) + + @property + def phase_components(self): + return self.component_list + + @property + def phase_funcs(self): + """List of all phase functions.""" + return self.get_quantity_funcs('phase_funcs_component') + + @property + def phase_deriv_funcs(self): + """List of derivative functions for phase components.""" + return self.get_deriv_funcs("deriv_funcs") + + @property + def d_phase_d_delay_funcs(self): + """List of d_phase_d_delay functions.""" + Dphase_Ddelay = [] + for cp in self.component_list: + Dphase_Ddelay += cp.phase_derivs_wrt_delay + return Dphase_Ddelay + + def phase(self, toas, delay=None, abs_phase=False): + """Return the model-predicted pulse phase for the given TOAs. + + Parameters + ---------- + toas: `~pint.toa.TOAs` object. + TOAs for evaluating the phase. + delay: `numpy.ndarray`, optional. + Input time delay values. If not given, phase will calculate the + delay. Default is None, + abs_phase: bool + Flag for using the absolute phase. Default is False. + + Return + ------ + `~pint.phase.Phase` object. The spin phase that at given TOAs. + """ + # First compute the delays to "pulsar time" + delay = self.delay(toas) + phase = Phase(np.zeros(toas.ntoas), np.zeros(toas.ntoas)) + # Then compute the relevant pulse phases + for pf in self.phase_funcs: + phase += Phase(pf(toas, delay)) + + # If the absolute phase flag is on, use the TZR parameters to compute + # the absolute phase. + if abs_phase: + if "AbsPhase" not in list(self.components.keys()): + # if no absolute phase (TZRMJD), add the component to the model and calculate it + from pint.models import absolute_phase + + self.add_component(absolute_phase.AbsPhase()) + self.make_TZR_toa( + toas + ) # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this? + tz_toa = self.get_TZR_toa(toas) + tz_delay = self.delay(tz_toa) + tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table))) + for pf in self.phase_funcs: + tz_phase += Phase(pf(tz_toa, tz_delay)) + return phase - tz_phase + else: + return phase + + def d_phase_d_param(self, toas, delay, param): + """Return the derivative of phase with respect to the parameter.""" + # TODO need to do correct chain rule stuff wrt delay derivs, etc + # Is it safe to assume that any param affecting delay only affects + # phase indirectly (and vice-versa)?? + par = getattr(self, param) + result = np.longdouble(np.zeros(toas.ntoas)) / par.units + param_phase_derivs = [] + phase_derivs = self.phase_deriv_funcs + if param in list(phase_derivs.keys()): + for df in phase_derivs[param]: + result += df(toas, param, delay).to( + result.unit, equivalencies=u.dimensionless_angles() + ) + else: + # Apply chain rule for the parameters in the delay. + # total_phase = Phase1(delay(param)) + Phase2(delay(param)) + # d_total_phase_d_param = d_Phase1/d_delay*d_delay/d_param + + # d_Phase2/d_delay*d_delay/d_param + # = (d_Phase1/d_delay + d_Phase2/d_delay) * + # d_delay_d_param + d_delay_d_p = self.d_delay_d_param(toas, param) + dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second + for dpddf in self.d_phase_d_delay_funcs: + dpdd_result += dpddf(toas, delay) + result = dpdd_result * d_delay_d_p + return result.to(result.unit, equivalencies=u.dimensionless_angles()) + + def d_phase_d_param_num(self, toas, param, step=1e-2): + """Return the derivative of phase with respect to the parameter. + + Compute the value numerically, using a symmetric finite difference. + + """ + # TODO : We need to know the range of parameter. + par = getattr(self, param) + ori_value = par.value + unit = par.units + if ori_value == 0: + h = 1.0 * step + else: + h = ori_value * step + parv = [par.value - h, par.value + h] + + phase_i = ( + np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled + ) + phase_f = ( + np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled + ) + for ii, val in enumerate(parv): + par.value = val + ph = self.phase(toas) + phase_i[:, ii] = ph.int + phase_f[:, ii] = ph.frac + res_i = -phase_i[:, 0] + phase_i[:, 1] + res_f = -phase_f[:, 0] + phase_f[:, 1] + result = (res_i + res_f) / (2.0 * h * unit) + # shift value back to the original value + par.quantity = ori_value + return result + + def d_phase_d_toa(self, toas, sample_step=None): + """Return the derivative of phase wrt TOA. + + Parameters + ---------- + toas : PINT TOAs class + The toas when the derivative of phase will be evaluated at. + sample_step : float optional + Finite difference steps. If not specified, it will take 1/10 of the + spin period. + + """ + copy_toas = copy.deepcopy(toas) + if sample_step is None: + pulse_period = 1.0 / (self.F0.quantity) + sample_step = pulse_period * 1000 + sample_dt = [-sample_step, 2 * sample_step] + + sample_phase = [] + for dt in sample_dt: + dt_array = [dt.value] * copy_toas.ntoas * dt._unit + deltaT = time.TimeDelta(dt_array) + copy_toas.adjust_TOAs(deltaT) + phase = self.phase(copy_toas) + sample_phase.append(phase) + # Use finite difference method. + # phase'(t) = (phase(t+h)-phase(t-h))/2+ 1/6*F2*h^2 + .. + # The error should be near 1/6*F2*h^2 + dp = sample_phase[1] - sample_phase[0] + d_phase_d_toa = dp.int / (2 * sample_step) + dp.frac / (2 * sample_step) + del copy_toas + return d_phase_d_toa.to(u.Hz) + + def get_spin_frequency(self, toas=None, apparent_frequency=False): + pass + + +class NoiseSector(ModelSector): + """ Class for holding all delay components and their APIs + """ + _methods = ('covariance_matrix_funcs', 'scaled_sigma_funcs', + 'basis_funcs', 'scaled_sigma', 'noise_model_designmatrix', + 'noise_model_basis_weight', 'noise_model_dimensions') + + def __init__(self, noise_components): + self.sector_name = 'NoiseComponent' + super(NoiseSector, self).__init__(noise_components) + + @property + def covariance_matrix_funcs(self,): + """List of covariance matrix functions.""" + return self.get_quantity_funcs('covariance_matrix_funcs') + + @property + def scaled_sigma_funcs(self,): + """List of scaled uncertainty functions.""" + return self.get_quantity_funcs('scaled_sigma_funcs') + + @property + def basis_funcs(self,): + """List of scaled uncertainty functions.""" + return self.get_quantity_funcs('basis_funcs') + + def noise_model_designmatrix(self, toas): + result = [] + if len(self.basis_funcs) == 0: + return None + + for nf in self.basis_funcs: + result.append(nf(toas)[0]) + return np.hstack([r for r in result]) + + def noise_model_basis_weight(self, toas): + result = [] + if len(self.basis_funcs) == 0: + return None + + for nf in self.basis_funcs: + result.append(nf(toas)[1]) + return np.hstack([r for r in result]) + + def noise_model_dimensions(self, toas): + """Returns a dictionary of correlated-noise components in the noise + model. Each entry contains a tuple (offset, size) where size is the + number of basis funtions for the component, and offset is their + starting location in the design matrix and weights vector.""" + result = {} + + # Correct results rely on this ordering being the + # same as what is done in the self.basis_funcs + # property. + if len(self.basis_funcs) > 0: + ntot = 0 + for nc in self.component_list: + bfs = nc.basis_funcs + if len(bfs) == 0: + continue + nbf = 0 + for bf in bfs: + nbf += len(bf(toas)[1]) + result[nc.category] = (ntot, nbf) + ntot += nbf + + return result + +builtin_sector_map = {'DelayComponent': DelaySector, + 'PhaseComponent': PhaseSector, + 'NoiseComponent': NoiseSector, + } From 4d36417d051409d2847ad983dbe62fe2fcbb4412 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Thu, 11 Jun 2020 19:01:46 -0400 Subject: [PATCH 10/16] black the code --- src/pint/models/model_sector.py | 90 ++++++++++++++++++++++----------- src/pint/models/timing_model.py | 63 +++++++++++++---------- 2 files changed, 97 insertions(+), 56 deletions(-) diff --git a/src/pint/models/model_sector.py b/src/pint/models/model_sector.py index 33b7389c4..c217e858f 100644 --- a/src/pint/models/model_sector.py +++ b/src/pint/models/model_sector.py @@ -24,23 +24,28 @@ class ModelSector(object): ---- The order of the component in the list is the order a component get computed. """ + _methods = tuple() def __init__(self, components): - if not hasattr(self, 'sector_name'): - raise TypeError("Please use ModelSector's subclass to " - "initialize for a specific model sector." - ) + if not hasattr(self, "sector_name"): + raise TypeError( + "Please use ModelSector's subclass to " + "initialize for a specific model sector." + ) # If only one component is given, convert it to a list if not isinstance(components, (list, tuple)): - components = [components,] + components = [ + components, + ] # Check if the components are the same type for cp in components: cp_type = get_component_type(cp) if cp_type != self.sector_name: - raise ValueError("Component {} is not a {} of" - " component.".format(cp.__class__.__name__, - self.sector_name)) + raise ValueError( + "Component {} is not a {} of" + " component.".format(cp.__class__.__name__, self.sector_name) + ) self.component_list = components self._parent = None @@ -92,12 +97,20 @@ def get_deriv_funcs(self, deriv_dict_name): class DelaySector(ModelSector): """ Class for holding all delay components and their APIs """ - _methods = ('delay_components','delay', 'delay_funcs', - 'get_barycentric_toas', 'd_delay_d_param', - 'delay_deriv_funcs', 'd_delay_d_param_num', 'delay') + + _methods = ( + "delay_components", + "delay", + "delay_funcs", + "get_barycentric_toas", + "d_delay_d_param", + "delay_deriv_funcs", + "d_delay_d_param_num", + "delay", + ) def __init__(self, delay_components, sector_map={}): - self.sector_name = 'DelayComponent' + self.sector_name = "DelayComponent" super(DelaySector, self).__init__(delay_components) @property @@ -106,7 +119,7 @@ def delay_components(self): @property def delay_funcs(self): - return self.get_quantity_funcs('delay_funcs_component') + return self.get_quantity_funcs("delay_funcs_component") @property def delay_deriv_funcs(self): @@ -174,7 +187,6 @@ def get_barycentric_toas(self, toas, cutoff_component=""): corr = self.delay(toas, cutoff_component, False) return tbl["tdbld"] * u.day - corr - def d_delay_d_param(self, toas, param, acc_delay=None): """Return the derivative of delay with respect to the parameter.""" par = getattr(self, param) @@ -229,12 +241,21 @@ class PhaseSector(ModelSector): Parameters ---------- """ - _methods = ('phase_components', 'phase_func', 'phase', 'phase_deriv_funcs', - 'd_phase_d_param', 'd_phase_d_param_num', 'd_phase_d_toa', - 'd_phase_d_delay_funcs', 'get_spin_frequency') + + _methods = ( + "phase_components", + "phase_func", + "phase", + "phase_deriv_funcs", + "d_phase_d_param", + "d_phase_d_param_num", + "d_phase_d_toa", + "d_phase_d_delay_funcs", + "get_spin_frequency", + ) def __init__(self, phase_components, sector_map={}): - self.sector_name = 'PhaseComponent' + self.sector_name = "PhaseComponent" super(PhaseSector, self).__init__(phase_components) @property @@ -244,7 +265,7 @@ def phase_components(self): @property def phase_funcs(self): """List of all phase functions.""" - return self.get_quantity_funcs('phase_funcs_component') + return self.get_quantity_funcs("phase_funcs_component") @property def phase_deriv_funcs(self): @@ -405,28 +426,35 @@ def get_spin_frequency(self, toas=None, apparent_frequency=False): class NoiseSector(ModelSector): """ Class for holding all delay components and their APIs """ - _methods = ('covariance_matrix_funcs', 'scaled_sigma_funcs', - 'basis_funcs', 'scaled_sigma', 'noise_model_designmatrix', - 'noise_model_basis_weight', 'noise_model_dimensions') + + _methods = ( + "covariance_matrix_funcs", + "scaled_sigma_funcs", + "basis_funcs", + "scaled_sigma", + "noise_model_designmatrix", + "noise_model_basis_weight", + "noise_model_dimensions", + ) def __init__(self, noise_components): - self.sector_name = 'NoiseComponent' + self.sector_name = "NoiseComponent" super(NoiseSector, self).__init__(noise_components) @property def covariance_matrix_funcs(self,): """List of covariance matrix functions.""" - return self.get_quantity_funcs('covariance_matrix_funcs') + return self.get_quantity_funcs("covariance_matrix_funcs") @property def scaled_sigma_funcs(self,): """List of scaled uncertainty functions.""" - return self.get_quantity_funcs('scaled_sigma_funcs') + return self.get_quantity_funcs("scaled_sigma_funcs") @property def basis_funcs(self,): """List of scaled uncertainty functions.""" - return self.get_quantity_funcs('basis_funcs') + return self.get_quantity_funcs("basis_funcs") def noise_model_designmatrix(self, toas): result = [] @@ -470,7 +498,9 @@ def noise_model_dimensions(self, toas): return result -builtin_sector_map = {'DelayComponent': DelaySector, - 'PhaseComponent': PhaseSector, - 'NoiseComponent': NoiseSector, - } + +builtin_sector_map = { + "DelayComponent": DelaySector, + "PhaseComponent": PhaseSector, + "NoiseComponent": NoiseSector, +} diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index c9b741843..fd62f6b97 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -18,9 +18,9 @@ from pint.models.parameter import strParameter, maskParameter from pint.phase import Phase from pint.utils import ( - PrefixError, - interesting_lines, - lines_of, + PrefixError, + interesting_lines, + lines_of, split_prefixed_name, get_component_type, ) @@ -140,7 +140,7 @@ class TimingModel(object): rather than to any particular component. """ - + def __init__(self, name="", components=[], external_sector_map={}): if not isinstance(name, str): raise ValueError( @@ -168,8 +168,9 @@ def __init__(self, name="", components=[], external_sector_map={}): ) for cp in components: - self.add_component(cp, validate=False, - external_sector_map=external_sector_map) + self.add_component( + cp, validate=False, external_sector_map=external_sector_map + ) def __repr__(self): return "{}(\n {}\n)".format( @@ -216,10 +217,12 @@ def __getattr__(self, name): errmsg = "'TimingModel' object and its component has no attribute" errmsg += " '%s'." % name if six.PY2: - sector_methods = super(TimingModel, self).__getattribute__('sector_methods') + sector_methods = super(TimingModel, self).__getattribute__( + "sector_methods" + ) else: - sector_methods = super().__getattribute__('sector_methods') - + sector_methods = super().__getattribute__("sector_methods") + if name in sector_methods.keys(): sector = sector_methods[name] return getattr(sector, name) @@ -244,9 +247,9 @@ def sector_methods(self): """ md = {} if six.PY2: - sectors = super(TimingModel, self).__getattribute__('model_sectors') + sectors = super(TimingModel, self).__getattribute__("model_sectors") else: - sectors = super().__getattribute__('model_sectors') + sectors = super().__getattribute__("model_sectors") for k, v in sectors.items(): for method in v._methods: md[method] = v @@ -373,30 +376,39 @@ def _make_sector(self, component, external_sector_map={}): all_sector_map = builtin_sector_map all_sector_map.update(external_sector_map) if not isinstance(component, (list, tuple)): - components = [component,] + components = [ + component, + ] # Assuem the first component's type is the type for all other components. com_type = get_component_type(components[0]) try: sector_cls = all_sector_map[com_type] except KeyError: - raise ValueError("Can not find the Sector class for" - " {}".format(com_type)) + raise ValueError("Can not find the Sector class for" " {}".format(com_type)) sector = sector_cls(component) if sector.sector_name in self.model_sectors.keys(): log.warn("Sector {} is already in the timing model. skip...") return common_method = set(sector._methods).intersection(self.sector_methods.keys()) if len(common_method) != 0: - raise ValueError("Sector {}'s methods have the same method with the current " - "sector methods. But sector methods should be unique." - "Please check .sector_methods for the current sector" - " methods.".format(sector.sector_name)) + raise ValueError( + "Sector {}'s methods have the same method with the current " + "sector methods. But sector methods should be unique." + "Please check .sector_methods for the current sector" + " methods.".format(sector.sector_name) + ) self.model_sectors[sector.sector_name] = sector # Link sector to the Timing model class. self.model_sectors[sector.sector_name]._parent = self - def add_component(self, component, order=DEFAULT_ORDER, force=False, - validate=True, external_sector_map={}): + def add_component( + self, + component, + order=DEFAULT_ORDER, + force=False, + validate=True, + external_sector_map={}, + ): """Add a component into TimingModel. Parameters @@ -432,8 +444,7 @@ def add_component(self, component, order=DEFAULT_ORDER, force=False, % component.__class__.__name__ ) else: - self._make_sector(component, - external_sector_map=external_sector_map) + self._make_sector(component, external_sector_map=external_sector_map) cur_cps = [] # link new component to TimingModel @@ -715,19 +726,19 @@ def covariance_matrix(self, toas): tbl = toas.table result = np.zeros((ntoa, ntoa)) # When there is no noise model. - if 'NoiseComponent'not in self.component_types: + if "NoiseComponent" not in self.component_types: result += np.diag(tbl["error"].quantity.to(u.s).value ** 2) return result for nf in self.covariance_matrix_funcs: result += nf(toas) return result - + @property def has_correlated_errors(self): """Whether or not this model has correlated errors.""" if "NoiseComponent" in self.component_types: - for nc in self.model_sectors['NoiseComponent'].component_list: + for nc in self.model_sectors["NoiseComponent"].component_list: # recursive if necessary if nc.introduces_correlated_errors: return True @@ -742,7 +753,7 @@ def scaled_sigma(self, toas): tbl = toas.table result = np.zeros(ntoa) * u.us # When there is no noise model. - if 'NoiseComponent' not in self.component_types: + if "NoiseComponent" not in self.component_types: result += tbl["error"].quantity return result From f8de9549ccbdf36bebb578d38f569508b82cf2dd Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Thu, 11 Jun 2020 19:14:10 -0400 Subject: [PATCH 11/16] update the example --- docs/examples/build_model_from_scratch.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/examples/build_model_from_scratch.md b/docs/examples/build_model_from_scratch.md index 69d149a9a..57db622ad 100644 --- a/docs/examples/build_model_from_scratch.md +++ b/docs/examples/build_model_from_scratch.md @@ -6,7 +6,7 @@ jupyter: extension: .md format_name: markdown format_version: '1.2' - jupytext_version: 1.4.1 + jupytext_version: 1.4.0 kernelspec: display_name: Python 3 language: python @@ -84,7 +84,7 @@ tm = TimingModel("NGC6400E", component_instances) To view all the components in `TimingModel` instance, we can use the property `.components`, which returns a dictionary (name as the key, component instance as the value). -Internally, the components are stored in a list(ordered list, you will see why this is important below) according to their types. All the delay type of components (subclasses of `DelayComponent` class) are stored in the `DelayComponent_list`, and the phase type of components (subclasses of `PhaseComponent` class) in the `PhaseComponent_list`. +Internally, the components are stored in the `.model_sectors`, a dictionary with the component type name as the key and `ModelSector` subclasses as the value. `ModelSector` class also have the methods to collect results from each component. For instance, `DelaySector` has `delay()` that compute the total delay from all the `DelayComponet` objects and `PhaseSector` has `phase()` for total phase. But these methods can be accessed from the `TimingModel` class. ```python # print the components in the timing model @@ -100,10 +100,10 @@ for (cp_name, cp_instance) in tm.components.items(): * `TimingModel.params()` : List all the parameters in the timing model from all components. * `TimingModel.setup()` : Setup the components (e.g., register the derivatives). * `TimingModel.validate()` : Validate the components see if the parameters are setup properly. -* `TimingModel.delay()` : Compute the total delay. -* `TimingModel.phase()` : Compute the total phase. -* `TimingModel.delay_funcs()` : List all the delay functions from all the delay components. -* `TimingModel.phase_funcs()` : List all the phase functions from all the phase components. +* `TimingModel.delay()` : Compute the total delay, if DelayComponent exists. +* `TimingModel.phase()` : Compute the total phase, if PhaseComponent exists. +* `TimingModel.delay_funcs()` : List all the delay functions from all the delay components, if DelayComponent exists. +* `TimingModel.phase_funcs()` : List all the phase functions from all the phase components, if DelayComponent exists. * `TimingModel.get_component_type()` : Get all the components from one category * `TimingModel.map_component()` : Get the component location. It returns the component's instance, order in the list, host list and its type. * `TimingModel.get_params_mapping()` : Report which component each parameter comes from. @@ -201,10 +201,10 @@ print("All components in timing model:") display(tm.components) print("\n") -print("Delay components in the DelayComponent_list (order matters!):") +print("Delay components in the DelayComponent sector (order matters!):") # print the delay component order, dispersion should be after the astrometry -display(tm.DelayComponent_list) +display(tm.model_sectors['DelayComponent'].componet_list) ``` The DM value can be set as we set the parameters above. From 7be596b39e5110328f66cf2232184fb027d1b5c2 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Fri, 12 Jun 2020 09:30:07 -0400 Subject: [PATCH 12/16] Fix bug of sector method name --- src/pint/models/model_sector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/model_sector.py b/src/pint/models/model_sector.py index c217e858f..1ec2fe492 100644 --- a/src/pint/models/model_sector.py +++ b/src/pint/models/model_sector.py @@ -244,7 +244,7 @@ class PhaseSector(ModelSector): _methods = ( "phase_components", - "phase_func", + "phase_funcs", "phase", "phase_deriv_funcs", "d_phase_d_param", From 801f812157d6ec5cf785408ca47b6992d498205c Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Fri, 12 Jun 2020 10:16:03 -0400 Subject: [PATCH 13/16] Fix a typo and bug --- docs/examples/build_model_from_scratch.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/examples/build_model_from_scratch.md b/docs/examples/build_model_from_scratch.md index 57db622ad..4a67fba64 100644 --- a/docs/examples/build_model_from_scratch.md +++ b/docs/examples/build_model_from_scratch.md @@ -84,7 +84,7 @@ tm = TimingModel("NGC6400E", component_instances) To view all the components in `TimingModel` instance, we can use the property `.components`, which returns a dictionary (name as the key, component instance as the value). -Internally, the components are stored in the `.model_sectors`, a dictionary with the component type name as the key and `ModelSector` subclasses as the value. `ModelSector` class also have the methods to collect results from each component. For instance, `DelaySector` has `delay()` that compute the total delay from all the `DelayComponet` objects and `PhaseSector` has `phase()` for total phase. But these methods can be accessed from the `TimingModel` class. +Internally, the components are stored in the `.model_sectors`, a dictionary with the component type name as the key and `ModelSector` subclasses as the value. `ModelSector` class also have the methods to collect results from each component. For instance, `DelaySector` has `delay()` that compute the total delay from all the `DelayComponent` objects and `PhaseSector` has `phase()` for total phase. But these methods can be accessed from the `TimingModel` class. ```python # print the components in the timing model @@ -204,7 +204,7 @@ print("\n") print("Delay components in the DelayComponent sector (order matters!):") # print the delay component order, dispersion should be after the astrometry -display(tm.model_sectors['DelayComponent'].componet_list) +display(tm.model_sectors['DelayComponent'].component_list) ``` The DM value can be set as we set the parameters above. From e226910b10a266926a9998fa3791294e6713dae8 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Fri, 12 Jun 2020 11:01:50 -0400 Subject: [PATCH 14/16] Fix bugs --- src/pint/models/model_sector.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/pint/models/model_sector.py b/src/pint/models/model_sector.py index 1ec2fe492..af7c3f89c 100644 --- a/src/pint/models/model_sector.py +++ b/src/pint/models/model_sector.py @@ -298,7 +298,11 @@ def phase(self, toas, delay=None, abs_phase=False): `~pint.phase.Phase` object. The spin phase that at given TOAs. """ # First compute the delays to "pulsar time" - delay = self.delay(toas) + if 'DelayComponent' in self.component_types: + delay = self.delay(toas) + else: + delay = np.zeros(toas.ntoas) * u.second + phase = Phase(np.zeros(toas.ntoas), np.zeros(toas.ntoas)) # Then compute the relevant pulse phases for pf in self.phase_funcs: From dd35e3a5dacf48f664216c708c5ba3144da47b34 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Fri, 12 Jun 2020 11:46:15 -0400 Subject: [PATCH 15/16] Fix more examples --- docs/examples/PINT_walkthrough.md | 2 +- docs/examples/understanding_timing_models.md | 18 +++++++++++------- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/docs/examples/PINT_walkthrough.md b/docs/examples/PINT_walkthrough.md index 50982e0e2..6a976c1f5 100644 --- a/docs/examples/PINT_walkthrough.md +++ b/docs/examples/PINT_walkthrough.md @@ -6,7 +6,7 @@ jupyter: extension: .md format_name: markdown format_version: '1.2' - jupytext_version: 1.4.1 + jupytext_version: 1.4.0 kernelspec: display_name: Python 3 language: python diff --git a/docs/examples/understanding_timing_models.md b/docs/examples/understanding_timing_models.md index 33109e4f4..ba81d56ab 100644 --- a/docs/examples/understanding_timing_models.md +++ b/docs/examples/understanding_timing_models.md @@ -6,7 +6,7 @@ jupyter: extension: .md format_name: markdown format_version: '1.2' - jupytext_version: 1.4.1 + jupytext_version: 1.4.0 kernelspec: display_name: Python 3 language: python @@ -56,12 +56,12 @@ The TimingModel class stores lists of the delay model components and phase compo ```python # When this list gets printed, it shows the parameters that are associated with each component as well. -m.DelayComponent_list +m.model_sectors['DelayComponent'].component_list ``` ```python # Now let's look at the phase components. These include the absolute phase, the spindown model, and phase jumps -m.PhaseComponent_list +m.model_sectors['PhaseComponent'].component_list ``` We can add a component to an existing model @@ -82,7 +82,7 @@ m.add_component(a, validate=False) ``` ```python -m.DelayComponent_list # The new instance is added to delay component list +m.model_sectors['DelayComponent'].component_list # The new instance is added to delay component list ``` There are two ways to remove a component from a model. This simplest is to use the `remove_component` method to remove it by name. @@ -112,7 +112,7 @@ from_list.remove(component) ``` ```python -m.DelayComponent_list # AstrometryEcliptic has been removed from delay list. +m.model_sectors['DelayComponent'].component_list # AstrometryEcliptic has been removed from delay list. ``` To switch the order of a component, just change the order of the component list @@ -121,7 +121,7 @@ To switch the order of a component, just change the order of the component list ```python # Let's look at the order of the components in the delay list first -_ = [print(dc.__class__) for dc in m.DelayComponent_list] +_ = [print(dc.__class__) for dc in m.model_sectors['DelayComponent'].component_list] ``` ```python @@ -133,7 +133,7 @@ from_list[order], from_list[new_order] = from_list[new_order], from_list[order] ```python # Print the classes to see the order switch -_ = [print(dc.__class__) for dc in m.DelayComponent_list] +_ = [print(dc.__class__) for dc in m.model_sectors['DelayComponent'].component_list] ``` Delays are always computed in the order of the DelayComponent_list @@ -181,3 +181,7 @@ for comp, tp in Component.component_types.items(): special ``` + +```python + +``` From fdf515fb0cab0b68cd98b433f7243e619062f5e5 Mon Sep 17 00:00:00 2001 From: Jing Luo Date: Fri, 12 Jun 2020 11:47:36 -0400 Subject: [PATCH 16/16] Black the code --- src/pint/fitter.py | 2 +- src/pint/models/model_sector.py | 2 +- src/pint/utils.py | 3 +-- tests/test_model_sector.py | 10 +++++----- 4 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index f5453f9de..7dcaad2c9 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -743,7 +743,7 @@ def fit_toas(self, maxiter=1, threshold=False, full_cov=False): accuracy where they both can be applied. """ chi2 = 0 - has_noise_model = 'NoiseComponent' in self.model.component_types + has_noise_model = "NoiseComponent" in self.model.component_types for i in range(max(maxiter, 1)): fitp = self.get_fitparams() fitpv = self.get_fitparams_num() diff --git a/src/pint/models/model_sector.py b/src/pint/models/model_sector.py index af7c3f89c..4e1e3d026 100644 --- a/src/pint/models/model_sector.py +++ b/src/pint/models/model_sector.py @@ -298,7 +298,7 @@ def phase(self, toas, delay=None, abs_phase=False): `~pint.phase.Phase` object. The spin phase that at given TOAs. """ # First compute the delays to "pulsar time" - if 'DelayComponent' in self.component_types: + if "DelayComponent" in self.component_types: delay = self.delay(toas) else: delay = np.zeros(toas.ntoas) * u.second diff --git a/src/pint/utils.py b/src/pint/utils.py index bef060af6..f3016f05b 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1417,8 +1417,7 @@ def get_component_type(component): comp_base = inspect.getmro(component.__class__) if comp_base[-2].__name__ != "Component": raise TypeError( - "Class '%s' is not a Component type class." - % component.__class__.__name__ + "Class '%s' is not a Component type class." % component.__class__.__name__ ) elif len(comp_base) < 3: raise TypeError( diff --git a/tests/test_model_sector.py b/tests/test_model_sector.py index d2c085d8e..ce2670536 100644 --- a/tests/test_model_sector.py +++ b/tests/test_model_sector.py @@ -23,15 +23,15 @@ def setup(self): self.all_components = Component.component_types def test_sector_init(self): - sector = DelaySector(self.all_components['AstrometryEquatorial']()) + sector = DelaySector(self.all_components["AstrometryEquatorial"]()) - assert sector.__class__.__name__ == 'DelaySector' + assert sector.__class__.__name__ == "DelaySector" assert len(sector.component_list) == 1 - assert hasattr(sector, 'delay') - assert hasattr(sector, 'delay_funcs') + assert hasattr(sector, "delay") + assert hasattr(sector, "delay_funcs") def test_copy(self): - sector = DelaySector(self.all_components['AstrometryEquatorial']()) + sector = DelaySector(self.all_components["AstrometryEquatorial"]()) copy_sector = deepcopy(sector) assert id(sector) != id(copy_sector)