diff --git a/setup.py b/setup.py index 7285e315d..a139bb00e 100644 --- a/setup.py +++ b/setup.py @@ -182,6 +182,7 @@ def EXTENSION(modname): 'xpsi.cellmesh.rays', 'xpsi.tools.energy_interpolator', 'xpsi.tools.energy_integrator', + 'xpsi.tools.energy_integrator_2Dedges', 'xpsi.tools.phase_integrator', 'xpsi.tools.phase_interpolator', 'xpsi.tools.synthesise', diff --git a/xpsi/Data.py b/xpsi/Data.py index da43d798f..1d7b8ea76 100644 --- a/xpsi/Data.py +++ b/xpsi/Data.py @@ -225,6 +225,8 @@ def trim_data( self , # Make the table of required channels assert min_channel >= self.channels[0] assert max_channel <= self.channels[-1] + if max_channel == -1: + max_channel = self.channels[-1] new_channels = [ min_channel <= c <= max_channel for c in self.channels] # Trim the counts and channels @@ -488,6 +490,8 @@ def from_pha( cls, path, counts_data = spectrum['RATE'] exposure = _np.double(1.0) + TLMIN= _np.int32(Header['TLMIN1']) + # No channels specified, use everything if channels is None: min_channel = _np.min( channel_data ) @@ -497,6 +501,20 @@ def from_pha( cls, path, min_channel = channels[0] max_channel = channels[-1] + if 'QUALITY' in spectrum.columns.names : + print('Bad bins detected...') + quality=spectrum['QUALITY'] + qualchannels = _np.where(quality==0) + spectrum = spectrum[qualchannels] + print(r'Spectrum reduced to {} channels'.format(_np.shape(qualchannels)[1])) + #channel_data = spectrum['CHANNEL'] + channel_data = _np.arange(TLMIN, TLMIN+len(spectrum['CHANNEL'])) + min_channel = _np.min( channel_data ) + max_channel = _np.max( channel_data ) + channels = _np.arange(min_channel,max_channel+1) + print(channel_data, channels, qualchannels) + counts_data = spectrum['COUNTS'] + # Get intrinsinc values first = 0 last = max_channel - min_channel @@ -622,4 +640,4 @@ def plot(self, num_rot = 2, dpi=200, colormap='inferno'): fig.colorbar( im , ax=ax3 , location='right', label='Counts') fig.set_dpi(dpi) - return fig, axs \ No newline at end of file + return fig, axs diff --git a/xpsi/Instrument.py b/xpsi/Instrument.py index a08894cb9..f4b6f30ee 100644 --- a/xpsi/Instrument.py +++ b/xpsi/Instrument.py @@ -79,6 +79,12 @@ class Instrument(ParameterSubspace): of edges must be :math:`p + 1`. The edges must be monotonically increasing. These edges will correspond to the nominal response matrix and any deviation from this matrix (see above). + + :param ndarray[p+1] quality_checks: + Check if spectra counts present good quality. Default is False. + + :param ndarray[p+1] goodchannels: + Instrument channel numbers for which we have good quality. :param tuple args: Container of parameter instances. @@ -89,14 +95,22 @@ class Instrument(ParameterSubspace): find its way to the base class. """ - def __init__(self, matrix, energy_edges, channels, channel_edges=None, + def __init__(self, matrix, energy_edges, channels, channel_edges=None, quality_checks=False, goodchannels=None, *args, **kwargs): self.matrix = matrix - self.energy_edges = energy_edges self.channels = channels + self._quality_checks = quality_checks + if self._quality_checks==True: + try: + #self._quality_checks = quality_checks + self._goodchannels = goodchannels + except TypeError: + raise ChannelError('If quality is checked, then a list of good channels must be returned') + self.energy_edges = energy_edges if channel_edges is not None: self.channel_edges = channel_edges + super(Instrument, self).__init__(*args, **kwargs) @@ -134,10 +148,13 @@ def matrix(self, matrix): try: for i in range(matrix.shape[0]): assert matrix[i,:].any() + except AssertionError: + raise ResponseError('Each row of the matrix must contain at least one positive number.') + try: for j in range(matrix.shape[1]): assert matrix[:,j].any() except AssertionError: - raise ResponseError('Each row and column of the matrix must contain at least one positive number.') + raise ResponseError('Each column of the matrix must contain at least one positive number.') self._matrix = matrix def construct_matrix(self): @@ -218,13 +235,27 @@ def energy_edges(self, energy_edges): except TypeError: raise EdgesError('Energy edges must be in a one-dimensional array of positive increasing values.') + try: - assert energy_edges.ndim == 1 - assert (energy_edges >= 0.0).all() - assert energy_edges.shape[0] == self._matrix.shape[1] + 1 - assert not (energy_edges[1:] <= energy_edges[:-1]).any() + #""" + if self._quality_checks==True: + if energy_edges.ndim == 2 : + for i in range(energy_edges.ndim): + assert (energy_edges[i] >= 0.0).all() + assert not (energy_edges[i, 1:] <= energy_edges[i, :-1]).any() + else: + assert energy_edges.ndim == 1 + assert (energy_edges >= 0.0).all() + assert energy_edges.shape[0] == self._matrix.shape[1] + 1 + assert not (energy_edges[1:] <= energy_edges[:-1]).any() + #""" + #assert energy_edges.ndim == 1 + #assert (energy_edges >= 0.0).all() + #assert energy_edges.shape[0] == self._matrix.shape[1] + 1 + #assert not (energy_edges[1:] <= energy_edges[:-1]).any() except AssertionError: - raise EdgesError('Energy edges must be in a one-dimensional array of positive increasing values, with a ' + raise EdgesError('Unless if there are bad channels causing the removal of response contribution from some energy bins,' + 'energy edges must be in a one-dimensional array of positive increasing values, with a ' 'length equal to number of energy intervals in the matrix + 1.') self._energy_edges = energy_edges @@ -266,12 +297,19 @@ def channel_edges(self, channel_edges): raise EdgesError('Channel edges must be in a one-dimensional array of positive increasing values.') try: - assert channel_edges.ndim == 1 - assert (channel_edges >= 0.0).all() - assert channel_edges.shape[0] == self._matrix.shape[0] + 1 - assert not (channel_edges[1:] <= channel_edges[:-1]).any() + if self._quality_checks==True: + assert channel_edges.ndim == 2 + for i in range(channel_edges.ndim): + assert (channel_edges[i] >= 0.0).all() + assert not (channel_edges[i, 1:] <= channel_edges[i, :-1]).any() + #assert channel_edges.shape[0] == self._matrix.shape[0] + 1 + else: + assert channel_edges.ndim == 1 + assert (channel_edges >= 0.0).all() + assert channel_edges.shape[0] == self._matrix.shape[0] + 1 + assert not (channel_edges[1:] <= channel_edges[:-1]).any() except AssertionError: - raise EdgesError('Channel edges must be in a one-dimensional array of positive increasing values, with a ' + raise EdgesError('Unless if there are bad channels, channel edges must be in a one-dimensional array of positive increasing values, with a ' 'length equal to the number of channel intervals in the matrix + 1.') self._channel_edges = channel_edges @@ -311,12 +349,39 @@ def channels(self, channel_array): yield + @property + def goodchannels(self): + """ Get the array of channels with good quality. + + """ + return self._goodchannels + + @goodchannels.setter + @make_verbose('Setting channels with good quality', + 'Good channels set') + def goodchannels(self, channel_array): + if not isinstance(channel_array, _np.ndarray): + try: + channel_array = _np.array(channel_array) + except TypeError: + raise ChannelError('Channel numbers must be in an array.') + + try: + assert channel_array.ndim == 1 + assert (channel_array >= 0).all() + except AssertionError: + raise ChannelError('Good channel numbers must be in a one-dimensional array of positive integers (including zero)') + + self._goodchannels = channel_array + + yield + @make_verbose('Trimming instrument response', 'Instrument response trimmed') def trim_response(self, min_channel=0, max_channel=-1, - threshold=1e-5 ): + threshold=1e-5): """ Trim the instrument response to the specified channel range. :param int min_channel: @@ -328,7 +393,7 @@ def trim_response(self, :param float threshold: The threshold value to use for trimming the instrument response. Channels / inputs with a total response below this value will be removed. - + """ # Make the table of required channels @@ -343,7 +408,7 @@ def trim_response(self, new_matrix = self.matrix[new_channels_indexes] empty_channels = _np.all( new_matrix <= threshold, axis=1) empty_inputs = _np.all( new_matrix <= threshold, axis=0) - + # Apply to matrix and channels directly self.matrix = new_matrix[~empty_channels][:,~empty_inputs] self.channels = self.channels[new_channels_indexes][ ~empty_channels ] @@ -352,8 +417,20 @@ def trim_response(self, new_energy_edges = [ self.energy_edges[k] for k in range(len(empty_inputs)) if not empty_inputs[k] ] self.energy_edges = _np.hstack( (new_energy_edges , self.energy_edges[ _np.where( self.energy_edges == new_energy_edges[-1] )[0] + 1 ] ) ) if hasattr( self , 'channel_edges' ): - new_channels_edges = [ self.channel_edges[ _np.where(old_channels==chan)[0][0]] for chan in self.channels] - self.channel_edges = _np.hstack( (new_channels_edges , self.channel_edges[_np.where( old_channels==self.channels[-1])[0] + 1]) ) + if hasattr(self, 'goodchannels'): + EMINS, EMAXS = self.channel_edges[0], self.channel_edges[1] + MINMAXS = _np.concatenate([EMINS,EMAXS]) + minmaxs = MINMAXS.tolist() + old_channels_edges = _np.sort(np.array(list(dict.fromkeys(minmaxs)))) + new_channels_edges = [ old.channel_edges[ _np.where(old_channels==chan)[0][0]] for chan in self.channels] + allchannel_edges = _np.hstack( (new_channels_edges , old.channel_edges[_np.where( old_channels==self.channels[-1])[0] + 1]) ) + EnewMINS = allchannel_edges[:-1][self.goodchannels] + EnewMAXS = allchannel_edges[1:][self.goodchannels] + self.channel_edges = _np.stack((EnewMINS,EnewMAXS)) + + else: + new_channels_edges = [ self.channel_edges[ _np.where(old_channels==chan)[0][0]] for chan in self.channels] + self.channel_edges = _np.hstack( (new_channels_edges , self.channel_edges[_np.where( old_channels==self.channels[-1])[0] + 1]) ) # Print if any trimming happens if empty_inputs.sum() > 0: @@ -379,6 +456,7 @@ def from_ogip_fits(cls, min_input=1, max_input=-1, datafolder=None, + Data_path=None, **kwargs): """ Loading method for Instrument using OGIP defined ARF/RMF or RSP. @@ -388,6 +466,9 @@ def from_ogip_fits(cls, :param str | None ARF_path: The path to the ARF file which should be OGIP compliant or None if the RMF_path points to a RSP file. + :param str | None Data_path: + The path to the Data file which should be OGIP compliant. Only to mention if there are quality information. + :param int min_channel: The minimum channel for which the instrument response is loaded. @@ -408,6 +489,7 @@ def from_ogip_fits(cls, """ if datafolder: + Data_path = _os.path.join( datafolder, Data_path ) if Data_path is not None else None ARF_path = _os.path.join( datafolder, ARF_path ) if ARF_path is not None else None RMF_path = _os.path.join( datafolder, RMF_path ) @@ -426,12 +508,12 @@ def from_ogip_fits(cls, assert RMF_instr == ARF_hdul['SPECRESP'].header['INSTRUME'] # Get the values and change the -1 values if requried + if min_channel == 0: + min_channel = TLMIN if max_channel == -1: max_channel = DETCHANS -1 if max_input == -1: max_input = NUMGRP - channels = _np.arange( min_channel, max_channel+1 ) - inputs = _np.arange( min_input, max_input+1 ) # Perform routine checks assert min_channel >= TLMIN and max_channel <= TLMAX @@ -442,6 +524,10 @@ def from_ogip_fits(cls, RMF_MATRIX = RMF_hdul['MATRIX'].data RMF_EBOUNDS = RMF_hdul['EBOUNDS'].data + # Get all the channels and input energies + channels = _np.array( RMF_EBOUNDS['CHANNEL'] ) + inputs = _np.arange( 1, NUMGRP+1 ) + # Get the channels from the data RMF = _np.zeros((DETCHANS, NUMGRP)) for i, (N_GRP, F_CHAN, N_CHAN, RMF_line) in enumerate( zip(RMF_MATRIX['N_GRP'], RMF_MATRIX['F_CHAN'], RMF_MATRIX['N_CHAN'], RMF_MATRIX['MATRIX']) ): @@ -462,18 +548,24 @@ def from_ogip_fits(cls, if n_chan == 0: continue - RMF[f_chan:f_chan+n_chan,i] += RMF_line[n_skip:n_skip+n_chan] + RMF[f_chan-TLMIN:f_chan+n_chan-TLMIN,i] += RMF_line[n_skip:n_skip+n_chan] n_skip += n_chan + # Get the indexes of channels + channels_indexes = _np.where( ( channels >= min_channel ) & ( channels <= max_channel ) )[0] + inputs_indexes = _np.where( ( inputs >= min_input ) & ( inputs <= max_input ) )[0] + channels = channels[channels_indexes] + inputs = inputs[inputs_indexes] + RMF = RMF[channels_indexes][:,inputs_indexes] + # Make the RSP, depending on the input files if ARF_path is None: - RSP = RMF[min_channel:max_channel+1,min_input-1:max_input] + RSP = RMF ARF_area = RMF.sum( axis=0 ) else: ARF = Table.read(ARF_path, 'SPECRESP') - ARF_area = ARF['SPECRESP'] + ARF_area = ARF['SPECRESP'][inputs_indexes] RSP = RMF * ARF_area - RSP = RSP[min_channel:max_channel+1,min_input-1:max_input] # Find empty columns and lines empty_channels = _np.all(RSP == 0, axis=1) @@ -490,18 +582,32 @@ def from_ogip_fits(cls, # Get the edges of energies for both input and channel energy_edges = _np.append( RMF_MATRIX['ENERG_LO'][inputs-1], RMF_MATRIX['ENERG_HI'][inputs[-1]-1]).astype(dtype=_np.double) - channel_energy_edges = _np.append(RMF_EBOUNDS['E_MIN'][channels],RMF_EBOUNDS['E_MAX'][channels[-1]]) + channel_energy_edges = _np.append(RMF_EBOUNDS['E_MIN'][channels-TLMIN],RMF_EBOUNDS['E_MAX'][channels[-1]-TLMIN]) + + # If required, check for quality. + quality_checks=False + goodchannels=None + if Data_path is not None: + with fits.open( Data_path ) as Data_hdul: + try: + QUAL = Data_hdul['SPECTRUM'].data['QUALITY'] + quality_checks=True + goodchannels=_np.where(QUAL==0) + except TypeError: + raise ChannelError('Quality of source channels should be specified!') # Make the instrument Instrument = cls(RSP, energy_edges, channels, channel_energy_edges, + quality_checks=quality_checks, + goodchannels=goodchannels, **kwargs) # Add ARF and RMF for plotting - Instrument.RMF = RMF[min_channel:max_channel+1,min_input-1:max_input][~empty_channels][:,~empty_inputs] - Instrument.ARF = ARF_area[min_input-1:max_input][~empty_inputs] + Instrument.RMF = RMF[~empty_channels][:,~empty_inputs] + Instrument.ARF = ARF_area[~empty_inputs] Instrument.name = RMF_instr return Instrument @@ -699,4 +805,4 @@ def from_ogip_fits(cls, pileup = XrayPileup(Instrument) Instrument.pileup = pileup - return Instrument \ No newline at end of file + return Instrument diff --git a/xpsi/PostProcessing/_corner.py b/xpsi/PostProcessing/_corner.py index 5b5f4876e..becc29560 100644 --- a/xpsi/PostProcessing/_corner.py +++ b/xpsi/PostProcessing/_corner.py @@ -499,7 +499,10 @@ def _plot_triangle(self, contour_args = self.get_attr('contours') contour_args.reverse() - if len(getdist_bcknds) == 1: + if "legend_labels" in kwargs: + legend_labels = kwargs.pop('legend_labels') + assert len(legend_labels) == len(getdist_bcknds), 'There must be as many legend labels as runs to plot.' + elif len(getdist_bcknds) == 1: legend_labels = None elif len(self._subset) > 1: legend_labels = self.get_attr('parent_ID') @@ -672,7 +675,8 @@ def _plot_triangle(self, bootstrap = bootstrap, n_simulate = kwargs.get('n_simulate'), force_draw = force_draw_i, - prior_samples_fname=prior_samples_fname) + prior_samples_fname=prior_samples_fname, + priors_identical=priors_identical) if (i==0 and priors_identical): break @@ -754,7 +758,8 @@ def _add_prior_density(self, plotter, posterior, KL_divergence, KL_base, bootstrap, n_simulate, force_draw, - prior_samples_fname): + prior_samples_fname, + priors_identical=False): """ Crudely estimate the prior density. Kullback-Leibler divergence estimated in bits for a combined run or @@ -790,6 +795,8 @@ def _add_prior_density(self, plotter, posterior, _np.save(samples_npy,samples) color, lw = (run.contours[key] for key in ('color', 'lw')) + if priors_identical: + color = 'black' quantiles = [None] * 3 diff --git a/xpsi/Signal.py b/xpsi/Signal.py index a74d96636..f97492ff4 100644 --- a/xpsi/Signal.py +++ b/xpsi/Signal.py @@ -8,6 +8,7 @@ from xpsi.Interstellar import Interstellar from xpsi.tools.energy_integrator import energy_integrator +from xpsi.tools.energy_integrator_2Dedges import energy_integrator_2Dedges #from xpsi.tools.energy_interpolator import energy_interpolator #from xpsi.tools.phase_integrator import phase_integrator @@ -17,6 +18,8 @@ from copy import deepcopy +from scipy import sparse + class LikelihoodError(xpsiError): """ Raised if there is a problem with the value of the log-likelihood. """ @@ -120,6 +123,19 @@ def __init__(self, else: self._instrument = instrument self._original_instrument = deepcopy( instrument ) + + if self._instrument._quality_checks == True: + #import pdb; pdb.set_trace() + newRSP=self._instrument.matrix[:, self._instrument.goodchannels[0]] + empty_inputs = _np.all( newRSP == 0, axis=1) + if (empty_inputs==True).any(): + newRSP = newRSP[:][~empty_inputs] + edgesmin=self._instrument.energy_edges[:-1][_np.where(empty_inputs==False)[0]] + edgesmax=self._instrument.energy_edges[1:][_np.where(empty_inputs==False)[0]] + self._instrument.energy_edges=_np.stack((edgesmin,edgesmax)) + self._instrument.matrix=newRSP + self._instrument.channels=self._data.channels + # Trimming the data and response so they fit together if min_channel != 0 or max_channel != -1: @@ -129,17 +145,26 @@ def __init__(self, print('WARNING : There is no counts in data object. This is normal if you are trying to synthesise data.' 'Otherwise something is very wrong and do not continue') self._instrument.trim_response(min_channel, max_channel) + if hasattr(self._instrument, 'pileup'): + self._instrument.pileup.instrument.trim_response(min_channel, max_channel) + self._instrument.pileup.arf_data = self._instrument.pileup.instrument.ARF + self._instrument.pileup.rmf_data = self._instrument.pileup.instrument.RMF + self._instrument.pileup.energies = self._instrument.pileup.instrument.energies + # Check that the channel arrays match a, b = data.index_range - if (len(self._data.channels) != len(self._instrument.channels[a:b])): - raise ChannelError( 'Size of the channel array declared for event data ' - 'does not match the size of the channel array declared' - ' for the loaded instrument response (sub)matrix.') + if ( len(self._data.channels) != (b-a) ): + raise ChannelError( 'Size of the channel array declared for event data does not match the declared size.') - if (self._data.channels != self._instrument.channels[a:b]).any(): - raise ChannelError('Channel array declared for event data does not ' - 'match channel array declared for the loaded ' - 'instrument response (sub)matrix.') + # Check that channels of instrument and data can be matched, + try: + a_instrument = _np.where( self._instrument.channels == self._data.channels[0] )[0][0] + b_instrument = _np.where( self._instrument.channels == self._data.channels[-1] )[0][0] + self._instrument_index_range_channels = ( a_instrument, b_instrument + 1 ) + assert not (self._data.channels != self._instrument.channels[a_instrument:b_instrument+1]).any() + except ChannelError or IndexError: + raise ChannelError('Channel array declared for event data does not match channel array declared for the loaded ' + 'instrument response (sub)matrix. The data channels need to be a subset of the instrument channels.') # Check that they come from the same instrument if hasattr( self._data , 'instrument' ) and hasattr( self._instrument , 'name' ): @@ -156,9 +181,14 @@ def __init__(self, self._background = None if support is not None: - if self._data.counts.shape[0]==support.shape[0]: self._support = support + elif self._instrument.channels.shape[0]==support.shape[0]: + self._support = support[a_instrument:b_instrument+1] + elif self._original_instrument.channels.shape[0]==support.shape[0]: + a_original_instrument = _np.where( self._original_instrument.channels == self._data.channels[0] )[0][0] + b_original_instrument = _np.where( self._original_instrument.channels == self._data.channels[-1] )[0][0] + self._support = support[a_original_instrument:b_original_instrument+1] else: raise TypeError("Data spectrum and background support must the have same shape") else: @@ -287,11 +317,6 @@ def instrument(self): def original_data(self): """ Get the a copy of the original instance of :class:`~.Data.Data`.""" return self._original_data - - @property - def original_instrument(self): - """ Get the a copy of the original instance of :class:`~.Instrument.Instrument`.""" - return self._original_instrument @property def photosphere(self): @@ -316,6 +341,7 @@ def _identify_waveband(self): """ a, b = self._data.index_range + # Now find the first non-zero energy inputs def search(i, j, k): while self._instrument.matrix[i,j] == 0.0: j += k @@ -327,6 +353,8 @@ def search(i, j, k): self._input_interval_range = (a, b) self._energy_edges = self._instrument.energy_edges[a:b + 1] self._energy_mids = (self._energy_edges[:-1] + self._energy_edges[1:])/2.0 + if self._energy_mids.ndim == 2: + self._energy_mids = self._energy_mids[0] @property def fast_energies(self): @@ -382,10 +410,16 @@ def register(self, signals, fast_mode=False, threads=1): if component is None: fast_total_counts.append(None) else: - integrated = energy_integrator(threads, - component, - _np.log10(self.fast_energies), - _np.log10(self._energy_edges)) + if self._energy_edges.ndim == 2: + integrated = energy_integrator_2Dedges(threads, + component, + _np.log10(self.fast_energies), + _np.log10(self._energy_edges)) + else: + integrated = energy_integrator(threads, + component, + _np.log10(self.fast_energies), + _np.log10(self._energy_edges)) # move interstellar to star? if self._interstellar is not None: @@ -393,7 +427,7 @@ def register(self, signals, fast_mode=False, threads=1): temp = self._instrument(integrated, self._input_interval_range, - self._data.index_range) + self._instrument_index_range_channels) fast_total_counts.append(_np.sum(temp)) @@ -432,11 +466,18 @@ def register(self, signals, fast_mode=False, threads=1): for hotRegion in signals: integrated = None + #import pdb; pdb.set_trace() for component in hotRegion: - temp = energy_integrator(threads, - component, - _np.log10(self._energies), - _np.log10(self._energy_edges)) + if self._energy_edges.ndim == 2: + temp = energy_integrator_2Dedges(threads, + component, + _np.log10(self._energies), + _np.log10(self._energy_edges)) + else: + temp = energy_integrator(threads, + component, + _np.log10(self._energies), + _np.log10(self._energy_edges)) try: integrated += temp except TypeError: @@ -445,12 +486,13 @@ def register(self, signals, fast_mode=False, threads=1): if self.cache: self.incident_flux_signals = integrated.copy() + #pdb.set_trace() if self._interstellar is not None: self._interstellar(self._energy_mids, integrated) self.signals = self._instrument(integrated, self._input_interval_range, - self._data.index_range) + self._instrument_index_range_channels) if self._background is not None: try: @@ -463,7 +505,7 @@ def register(self, signals, fast_mode=False, threads=1): self._background.registered_background = \ self._instrument(self._background.incident_background, self._input_interval_range, - self._data.index_range) + self._instrument_index_range_channels) @property def num_components(self): @@ -782,10 +824,17 @@ def construct_energy_array(num_energies, signals, max_energy=None): try: MAX except NameError: - MAX = signal.energy_edges[-1] + if signal.energy_edges.ndim==2: + MAX = signal.energy_edges[1][-1] + else: + MAX = signal.energy_edges[-1] s = signal - E = signal.energy_edges[-1] + if signal.energy_edges.ndim==2: + E = signal.energy_edges[1][-1] + else: + E = signal.energy_edges[-1] + if E > MAX: MAX = E s = signal @@ -798,9 +847,16 @@ def construct_energy_array(num_energies, signals, max_energy=None): try: MIN except NameError: - MIN = signal.energy_edges[0] + if signal.energy_edges.ndim==2: + MIN = signal.energy_edges[0][0] + else: + MIN = signal.energy_edges[0] - E = signal.energy_edges[0] + if signal.energy_edges.ndim==2: + E = signal.energy_edges[0][0] + else: + E = signal.energy_edges[0] + if E < MIN: MIN = E @@ -810,7 +866,11 @@ def construct_energy_array(num_energies, signals, max_energy=None): del MAX # find global limits - _signal_max = ordered[0].energy_edges[-1] + if ordered[0].energy_edges.ndim==2: + _signal_max = ordered[0].energy_edges[1][-1] + else: + _signal_max = ordered[0].energy_edges[-1] + if max_energy is not None and max_energy < _signal_max: MAX = max_energy @@ -832,9 +892,16 @@ def construct_energy_array(num_energies, signals, max_energy=None): try: MIN except NameError: - MIN = signal.energy_edges[0] + if signal.energy_edges.ndim==2: + MIN = signal.energy_edges[0][0] + else: + MIN = signal.energy_edges[0] - E = signal.energy_edges[0] + if signal.energy_edges.ndim==2: + E = signal.energy_edges[0][0] + else: + E = signal.energy_edges[0] + if E < MIN: MIN = E diff --git a/xpsi/cellmesh/integrator.pxd b/xpsi/cellmesh/integrator.pxd index f23556b8e..c36047fd6 100644 --- a/xpsi/cellmesh/integrator.pxd +++ b/xpsi/cellmesh/integrator.pxd @@ -4,6 +4,7 @@ from GSL cimport (gsl_interp_eval, gsl_interp_accel, gsl_interp_accel_alloc, gsl_interp_steffen, + gsl_interp_linear, gsl_interp, gsl_interp_init, gsl_interp_free, diff --git a/xpsi/tools/core.pxd b/xpsi/tools/core.pxd index 4bff3de24..91b71ea94 100644 --- a/xpsi/tools/core.pxd +++ b/xpsi/tools/core.pxd @@ -2,6 +2,7 @@ from GSL cimport (gsl_interp_type, gsl_interp_akima_periodic, gsl_interp_akima, gsl_interp_steffen, + gsl_interp_linear, gsl_interp_cspline_periodic, gsl_interp_cspline) diff --git a/xpsi/tools/core.pyx b/xpsi/tools/core.pyx index e4e785717..df4b8205a 100644 --- a/xpsi/tools/core.pyx +++ b/xpsi/tools/core.pyx @@ -11,6 +11,7 @@ from .phase_integrator import phase_integrator from .phase_interpolator import phase_interpolator from .energy_integrator import energy_integrator +from .energy_integrator_2Dedges import energy_integrator_2Dedges from .energy_interpolator import energy_interpolator from .synthesise import synthesise_exposure @@ -18,7 +19,7 @@ from .synthesise import synthesise_given_total_count_number from libc.math cimport fabs -__interpolants__ = {'Akima': 0, 'Steffen': 1, 'Cubic': 2} +__interpolants__ = {'Akima': 0, 'Steffen': 1, 'Cubic': 2, 'Linear': 3} def get_phase_interpolant(): """ Get the globally-set phase interpolant. """ @@ -36,7 +37,7 @@ def set_phase_interpolant(interpolant): :param str interpolant: Name of the spline interpolant. Options are "Akima" (default), "Steffen" - (pre-v0.6 choice), and "Cubic". The first and last have periodic + (pre-v0.6 choice), "Linear" and "Cubic". The first and last have periodic boundary conditions. """ @@ -67,7 +68,7 @@ def set_energy_interpolant(interpolant): :param str interpolant: Name of the spline interpolant. Options are "Akima" (default), "Steffen" - (pre-v0.6 choice), and "Cubic". All have natural boundary conditions. + (pre-v0.6 choice), "Linear" and "Cubic". All have natural boundary conditions. """ global __energy_interpolant__ @@ -95,6 +96,8 @@ cdef const gsl_interp_type* _get_phase_interpolant() except *: return gsl_interp_steffen elif __phase_interpolant__ == 2: return gsl_interp_cspline_periodic + elif __phase_interpolant__ == 3: + return gsl_interp_linear else: raise ValueError('Invalid phase interpolant setting.') @@ -112,8 +115,10 @@ cdef const gsl_interp_type* _get_energy_interpolant() except *: return gsl_interp_steffen elif __energy_interpolant__ == 2: return gsl_interp_cspline + elif __energy_interpolant__ == 3: + return gsl_interp_linear else: - raise ValueError('Invalid eneergy interpolant setting.') + raise ValueError('Invalid energy interpolant setting.') cdef bint are_equal(double x, double y, double epsilon = 1.0e-12) noexcept nogil: if(fabs(x - y) < epsilon): diff --git a/xpsi/tools/energy_integrator_2Dedges.pyx b/xpsi/tools/energy_integrator_2Dedges.pyx new file mode 100644 index 000000000..fbdbaff3b --- /dev/null +++ b/xpsi/tools/energy_integrator_2Dedges.pyx @@ -0,0 +1,114 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +#cython: embedsignature=True + +import numpy as np +from cython.parallel cimport * +from libc.stdlib cimport malloc, free +from libc.math cimport pow, log + +from GSL cimport (gsl_interp, + gsl_interp_alloc, + gsl_interp_init, + gsl_interp_free, + gsl_interp_eval, + gsl_interp_eval_integ, + gsl_interp_accel, + gsl_interp_accel_alloc, + gsl_interp_accel_free, + gsl_interp_accel_reset) + +ctypedef gsl_interp_accel accel + +from .core cimport _get_phase_interpolant, gsl_interp_type + +def energy_integrator_2Dedges(size_t N_Ts, + double[:,::1] signal, + double[::1] energies, + double[:,::1] energy_edges): + """ Integrate a signal over energy intervals. + + :param size_t N_Ts: + Number of OpenMP threads to spawn. + + :param double[:,::1] signal: + A C-contiguous :class:`numpy.ndarray` of an energy-resolved (specific + flux) signal. Energy increases with row number. + + :param double[::1] energies: + A :class:`numpy.ndarray` of the logarithms (base 10) of the energies. + + :param double[:,::1] energy_edges: + A :class:`numpy.ndarray` of the logarithm (base 10) of the energy + interval edges. + + :returns: + A 2D :class:`numpy.ndarray` of the signal integrated over energy + intervals. Energy interval number increases with row number. + + """ + cdef const gsl_interp_type *_interpolant + + _interpolant = _get_phase_interpolant() + + cdef: + signed int ii + size_t i, j, T + double *cpy + double upper_energy + double max_energy = energies[energies.shape[0] - 1] + + double **_signal = malloc(sizeof(double*) * N_Ts) + gsl_interp **interp = malloc(sizeof(gsl_interp*) * N_Ts) + accel **acc = malloc(N_Ts * sizeof(accel*)) + + double[:,::1] binned_signal = np.zeros((signal.shape[1], + energy_edges.shape[1]), + dtype = np.double) + + for T in range(N_Ts): + acc[T] = gsl_interp_accel_alloc() + interp[T] = gsl_interp_alloc(_interpolant, energies.shape[0]) + _signal[T] = malloc(sizeof(double) * energies.shape[0]) + + for ii in prange(signal.shape[1], + nogil = True, + schedule = 'static', + num_threads = N_Ts, + chunksize = 1): + i = ii + T = threadid() + cpy = _signal[T] + + for j in range(energies.shape[0]): + cpy[j] = pow(10.0, energies[j]) * signal[j,i] * log(10.0) + + gsl_interp_accel_reset(acc[T]) + gsl_interp_init(interp[T], &(energies[0]), cpy, energies.shape[0]) + + for j in range(energy_edges.shape[1] ): + if energy_edges[1][j] > max_energy: + upper_energy = max_energy + else: + upper_energy = energy_edges[1][j] + binned_signal[i,j] = gsl_interp_eval_integ(interp[T], + &(energies[0]), + cpy, + energy_edges[0][j], + upper_energy, + acc[T]) + if energy_edges[1][j] > max_energy: + break + + for T in range(N_Ts): + gsl_interp_accel_free(acc[T]) + gsl_interp_free(interp[T]) + free(_signal[T]) + + free(_signal) + free(interp) + free(acc) + + return np.asarray(binned_signal.T, dtype = np.double, order = 'C') diff --git a/xpsi/utilities/BackgroundTools.py b/xpsi/utilities/BackgroundTools.py index 256a1524f..cb52cc1c8 100644 --- a/xpsi/utilities/BackgroundTools.py +++ b/xpsi/utilities/BackgroundTools.py @@ -181,8 +181,12 @@ def plotBackgroundSpectrum( XPSI_model, print( f"Maximum counts in an energy bin : {np.max( HotRegion_spectra.sum(axis=0) )}") # Extract channels - x0 = signal.instrument.channel_edges - x0 = ( x0[:-1] + x0[1:] ) / 2 + if signal.instrument.channel_edges.ndim == 2: + x01, x02 = signal.instrument.channel_edges[0], signal.instrument.channel_edges[1] + x0 = (x01+x02)/2 + else: + x0 = signal.instrument.channel_edges + x0 = ( x0[:-1] + x0[1:] ) / 2 # Extract background from samples if plot_range: @@ -225,8 +229,8 @@ def plotBackgroundSpectrum( XPSI_model, ax.plot(x0,Data_Spectrum,'--', color=mycolors[5], lw=2, label='Data light curve') # Plotting background support - if plot_support and signal.support is not None and signal.support[signal.support>0].any(): - support = signal.support * data.exposure_time + if plot_support and signal._support is not None and signal._support[signal._support>0].any(): + support = signal._support * data.exposure_time ax.fill_between(x0, support[:,0], support[:,1], color='red', alpha = 0.2, label='BKG prior support') # Finish the plot