From aeb3ac4769e1544e4e854c01feaf268bd52d266e Mon Sep 17 00:00:00 2001 From: yeganer Date: Wed, 21 Jun 2017 15:41:18 +0200 Subject: [PATCH 1/3] [TEP005][QOL] Add MontecarloRunner.spectrum_integrated This is a shortcut for obtaining an integrated spectrum for the wavelength grid defined in the configuration file. If incompatible configuration settings are detected, an IntegrationError is raised. --- scripts/tardis | 27 +++++------------ tardis/montecarlo/base.py | 7 +++++ tardis/montecarlo/formal_integral.py | 45 +++++++++++++++++++++++++++- 3 files changed, 58 insertions(+), 21 deletions(-) diff --git a/scripts/tardis b/scripts/tardis index e990c8a6671..1797aaa9669 100755 --- a/scripts/tardis +++ b/scripts/tardis @@ -67,28 +67,13 @@ tardis_config = config_reader.Configuration.from_yaml(args.config_fname) simulation = Simulation.from_config(tardis_config) -def get_integrated_spectrum(): - if ( - tardis_config.plasma.disable_electron_scattering and - tardis_config.plasma.line_interaction_type == 'downbranch' - ): - return simulation.runner.integrator.calculate_spectrum( - simulation.runner.spectrum_frequency[:-1], - 1000) - else: - raise config_reader.ConfigurationError( - 'The formal integral module currently works only for ' - 'plasma settings "downbranch" and disabled electron scattering' - '.\nWe got: {}'.format(str(tardis_config.plasma))) - - def get_virtual_spectrum(): # Catch warning when acessing invalid spectrum_virtual with warnings.catch_warnings(record=True) as w: spectrum = simulation.runner.spectrum_virtual if len(w) > 0 and w[-1]._category_name == 'UserWarning': - print( + warnings.warn( 'Virtual spectrum is not available, using the ' 'real packet spectrum instead.') spectrum = simulation.runner.spectrum @@ -99,15 +84,17 @@ print('Saving the {} spectrum.'.format(tardis_config.spectrum.method)) # Believe it or not, that is a fancy switch-case statement # We need functions so that only one spectrum is accessed so we can catch # the warning properly -spectrum = { +get_spectrum = { 'real': lambda: simulation.runner.spectrum, 'virtual': get_virtual_spectrum, - 'integrated': get_integrated_spectrum, + 'integrated': lambda: simulation.runner.spectrum_integrated, }[tardis_config.spectrum.method] if args.gdb: - import os; print(os.getpid()); raw_input() # Workaround to attach gdb + import os + print(os.getpid()) + raw_input() # Workaround to attach gdb if args.profile: import cProfile cProfile.runctx('simulation.run()', locals(), @@ -115,4 +102,4 @@ if args.profile: else: simulation.run() -spectrum().to_ascii(args.spectrum_fname) +get_spectrum().to_ascii(args.spectrum_fname) diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index 5c4ac7cb005..3b5df95db45 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -124,6 +124,13 @@ def spectrum_virtual(self): self.spectrum_frequency, self.montecarlo_virtual_luminosity) + @property + def spectrum_integrated(self): + return self.integrator.calculate_spectrum( + self.spectrum_frequency[:-1], + 1000, + raises=True) + @property def integrator(self): if self._integrator is None: diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 961344576c3..06b9f156b59 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -1,3 +1,4 @@ +import warnings import numpy as np import pandas as pd from astropy import units as u @@ -6,6 +7,9 @@ from tardis.montecarlo.spectrum import TARDISSpectrum +class IntegrationError(Exception): + pass + class FormalIntegrator(object): def __init__(self, model, plasma, runner): @@ -13,10 +17,49 @@ def __init__(self, model, plasma, runner): self.plasma = plasma self.runner = runner - def calculate_spectrum(self, frequency, N=1000): + def check(self, raises=False): + ''' + A method that determines if the formal integral can be performed with + the current configuration settings + + The function returns False if the configuration conflicts with the + required settings. If raises evaluates to True, then a + IntegrationError is raised instead + ''' + def raise_or_return(message): + if raises: + raise IntegrationError(message) + else: + warnings.warn(message) + return False + + for obj in (self.model, self.plasma, self.runner): + if obj is None: + return raise_or_return( + 'The integrator is missing either model, plasma or ' + 'runner. Please make sure these are provided to the ' + 'FormalIntegrator.' + ) + + if not self.runner.line_interaction_type == 'downbranch': + return raise_or_return( + 'The FormalIntegrator currently only works for ' + 'line_interaction_type == "downbranch"' + ) + + if self.runner.sigma_thomson.value > 1e-100: + return raise_or_return( + 'The FormalIntegrator currently only works with ' + 'disable_electron_scattering turned on.' + ) + + return True + + def calculate_spectrum(self, frequency, N=1000, raises=False): # Very crude implementation # The c extension needs bin centers (or something similar) # while TARDISSpectrum needs bin edges + self.check(raises) frequency = frequency.to('Hz', u.spectral()) luminosity = u.Quantity( formal_integral( From 22e15c5448de0e1f98568c07604d2827efc2688c Mon Sep 17 00:00:00 2001 From: yeganer Date: Thu, 22 Jun 2017 09:13:06 +0200 Subject: [PATCH 2/3] [TEP005] Change Integrator to raise by default --- tardis/montecarlo/base.py | 3 +-- tardis/montecarlo/formal_integral.py | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index 3b5df95db45..eca0abdba1f 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -128,8 +128,7 @@ def spectrum_virtual(self): def spectrum_integrated(self): return self.integrator.calculate_spectrum( self.spectrum_frequency[:-1], - 1000, - raises=True) + 1000) @property def integrator(self): diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 06b9f156b59..b34d213c7d8 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -17,7 +17,7 @@ def __init__(self, model, plasma, runner): self.plasma = plasma self.runner = runner - def check(self, raises=False): + def check(self, raises=True): ''' A method that determines if the formal integral can be performed with the current configuration settings @@ -55,7 +55,7 @@ def raise_or_return(message): return True - def calculate_spectrum(self, frequency, N=1000, raises=False): + def calculate_spectrum(self, frequency, N=1000, raises=True): # Very crude implementation # The c extension needs bin centers (or something similar) # while TARDISSpectrum needs bin edges From 075951e1bce1436b48a154c902516edbfe823bb7 Mon Sep 17 00:00:00 2001 From: yeganer Date: Thu, 22 Jun 2017 13:46:31 +0200 Subject: [PATCH 3/3] Add default value for number of integration points --- tardis/montecarlo/base.py | 3 +-- tardis/montecarlo/formal_integral.py | 10 ++++++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index eca0abdba1f..0a618ab75b4 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -127,8 +127,7 @@ def spectrum_virtual(self): @property def spectrum_integrated(self): return self.integrator.calculate_spectrum( - self.spectrum_frequency[:-1], - 1000) + self.spectrum_frequency[:-1]) @property def integrator(self): diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index b34d213c7d8..b9200500c28 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -10,12 +10,14 @@ class IntegrationError(Exception): pass + class FormalIntegrator(object): - def __init__(self, model, plasma, runner): + def __init__(self, model, plasma, runner, points=1000): self.model = model self.plasma = plasma self.runner = runner + self.points = points def check(self, raises=True): ''' @@ -55,12 +57,14 @@ def raise_or_return(message): return True - def calculate_spectrum(self, frequency, N=1000, raises=True): + def calculate_spectrum(self, frequency, points=None, raises=True): # Very crude implementation # The c extension needs bin centers (or something similar) # while TARDISSpectrum needs bin edges self.check(raises) + N = points or self.points frequency = frequency.to('Hz', u.spectral()) + luminosity = u.Quantity( formal_integral( self, @@ -68,6 +72,7 @@ def calculate_spectrum(self, frequency, N=1000, raises=True): N), 'erg' ) * (frequency[1] - frequency[0]) + # Ugly hack to convert to 'bin edges' frequency = u.Quantity( np.concatenate([ @@ -76,6 +81,7 @@ def calculate_spectrum(self, frequency, N=1000, raises=True): frequency.value[-1] + np.diff(frequency.value)[-1] ]]), frequency.unit) + return TARDISSpectrum( frequency, luminosity