From 74e8efd3beeb687de1dac19ec01c78451487a317 Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 11 Oct 2023 09:23:07 -0700 Subject: [PATCH 1/8] better reporting of tracebacks --- fitburst/analysis/fitter.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/fitburst/analysis/fitter.py b/fitburst/analysis/fitter.py index d276521..b0e9ee8 100644 --- a/fitburst/analysis/fitter.py +++ b/fitburst/analysis/fitter.py @@ -10,6 +10,7 @@ from scipy.optimize import least_squares import fitburst.routines.derivative as deriv import numpy as np +import traceback import sys class LSFitter: @@ -297,7 +298,7 @@ def fit(self, exact_jacobian: bool = True) -> None: except Exception as exc: print("ERROR: solver encountered a failure! Debug!") - print(sys.exc_info()) + print(traceback.format_exc()) def fix_parameter(self, parameter_list: list) -> None: """ @@ -474,6 +475,7 @@ def _compute_fit_statistics(self, spectrum_observed: float, fit_result: object) except Exception as exc: print(f"ERROR: {exc}; designating fit as unsuccessful...") + print(traceback.format_exc()) self.success = False def _set_weights(self) -> None: From 06ee7d5960d843b1f856b04cfee994237e3fa203 Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 11 Oct 2023 09:39:12 -0700 Subject: [PATCH 2/8] redefined PBF to drop the sc-time parameter and factor of 0.5 from the amplitude. --- fitburst/routines/profile.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/fitburst/routines/profile.py b/fitburst/routines/profile.py index 889ae05..e190eff 100644 --- a/fitburst/routines/profile.py +++ b/fitburst/routines/profile.py @@ -48,8 +48,8 @@ def compute_profile_gaussian(values: float, mean: float, width: float, return profile -def compute_profile_pbf(times: np.ndarray, toa: float, width: float, - sc_time: float) -> float: +def compute_profile_pbf(time: float, toa: float, width: float, freq: float, ref_freq: float, + sc_time_ref: float, sc_index: float = -4.) -> float: """ Computes a one-dimensional pulse broadening function (PBF) using the analytical solution of a Gaussian profile convolved with a one-side @@ -57,8 +57,8 @@ def compute_profile_pbf(times: np.ndarray, toa: float, width: float, Parameters ---------- - times : array_like - an array of times at which to evaluate PBF + time : array_like + a value or array of times at which to evaluate PBF toa : float the time of arrival of the burst @@ -66,8 +66,17 @@ def compute_profile_pbf(times: np.ndarray, toa: float, width: float, width : float the temporal width of the burst - sc_time : float - the scattering timescale quantifying the one-sided exponential tail + freq : float + the electromagnetic frequency at which to evaluate the PBF + + ref_freq : float + the electromagnetic frequency at which to reference the PBF + + sc_time_ref : float + the scattering timescale at the frequency 'ref_freq' + + sc_index : float, optional + the exponent of the frequency dependence for the PBF Returns ------- @@ -78,10 +87,14 @@ def compute_profile_pbf(times: np.ndarray, toa: float, width: float, The PBF equation is taken from McKinnon et al. (2004). """ - amp_term = 1. / 2. / sc_time - exp_term_1 = np.exp(width**2 / 2 / sc_time**2) - exp_term_2 = np.exp(-(times - toa) / sc_time[:, None]) - erf_term = 1 + ss.erf((times - (toa + width**2 / sc_time[:, None])) / width / np.sqrt(2)) - profile = amp_term[:, None] * exp_term_1[:, None] * exp_term_2 * erf_term + # evaluate the separate terms that define the PBF. + amp_term = (freq / ref_freq) ** (-sc_index) + sc_time = sc_time_ref * (freq / ref_freq) ** sc_index + arg_exp = width ** 2 / 2 / sc_time ** 2 - (time - toa) / sc_time + exp_term = np.exp(arg_exp) + erf_term = 1 + ss.erf((time - (toa + width ** 2 / sc_time)) / width / np.sqrt(2)) + + # finally, compute the PBF. + profile = amp_term * exp_term * erf_term return profile From afa0994567ac7f6adbd6db8912a10d9f734ed39d Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 11 Oct 2023 09:43:05 -0700 Subject: [PATCH 3/8] rejigged instantiation of PBF function call --- fitburst/analysis/model.py | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/fitburst/analysis/model.py b/fitburst/analysis/model.py index 3db87fc..38a226d 100644 --- a/fitburst/analysis/model.py +++ b/fitburst/analysis/model.py @@ -216,8 +216,11 @@ def compute_model(self, data: float = None) -> float: current_profile = self.compute_profile( current_times_arr, 0.0, # since 'current_times' is already corrected for DM. - current_sc_time_scaled, + current_sc_time, + current_sc_idx, current_width, + current_freq_arr, + current_ref_freq, is_folded = self.is_folded, ) @@ -258,8 +261,8 @@ def compute_model(self, data: float = None) -> float: return np.sum(self.spectrum_per_component, axis=2) - def compute_profile(self, times: float, arrival_time: float, sc_time: float, - width: float, is_folded: bool = False) -> float: + def compute_profile(self, times: float, arrival_time: float, sc_time_ref: float, sc_index: float, + width: float, freqs: float, ref_freq: float, is_folded: bool = False) -> float: """ Returns the temporal profile, depending on input values of width and scattering timescale. @@ -272,12 +275,21 @@ def compute_profile(self, times: float, arrival_time: float, sc_time: float, arrival_time : float The arrival time of the burst - sc_time : float + sc_time_ref : float The scattering timescale of the burst (which depends on frequency label) + sc_index : float + The index of frequency dependence on the scattering timescale + width : float The intrinsic temporal width of the burst + freqs : float + The index of frequency dependence on the scattering timescale + + ref_freq : float + The index of frequency dependence on the scattering timescale + is_folded : bool, optional If true, then the temporal profile is computed over two realizations and then averaged down to one (in order to allow for wrapping of a folded pulse shape) @@ -304,6 +316,7 @@ def compute_profile(self, times: float, arrival_time: float, sc_time: float, # compute either Gaussian or pulse-broadening function, depending on inputs. profile = np.zeros(times_copy.shape, dtype=float) + sc_time = sc_time_ref * (freqs / ref_freq) ** sc_index if np.any(sc_time < np.fabs(0.15 * width)): profile = rt.profile.compute_profile_gaussian(times_copy, arrival_time, width) @@ -313,7 +326,9 @@ def compute_profile(self, times: float, arrival_time: float, sc_time: float, # floating-point overlow in the exp((-times - toa) / sc_time) term in the # PBF call. TODO: use a better, more transparent method for avoiding this. times_copy[times_copy < -5 * width] = -5 * width - profile = rt.profile.compute_profile_pbf(times_copy, arrival_time, width, sc_time) + profile = rt.profile.compute_profile_pbf( + times_copy, arrival_time, width, freqs, ref_freq, sc_time_ref, sc_index=sc_index + ) # if data are folded and time/profile data contain two realizations, then # average along the appropriate axis to obtain a single realization. From 20b1276cf481d5d0bc0f5550103e18b9f61b1db0 Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 11 Oct 2023 09:44:32 -0700 Subject: [PATCH 4/8] first round of overhaul for derivative code --- fitburst/routines/derivative.py | 959 +++++++++++++++++++++----------- 1 file changed, 643 insertions(+), 316 deletions(-) diff --git a/fitburst/routines/derivative.py b/fitburst/routines/derivative.py index 6e1e86f..e14a293 100644 --- a/fitburst/routines/derivative.py +++ b/fitburst/routines/derivative.py @@ -7,9 +7,289 @@ fit parameters. """ +from ..backend import general import scipy.special as ss import numpy as np +def argument_erf(burst_width: float, time_diff: float, scattering_timescale: float) -> float: + """ + Computes the argument of the error function in the scatter-broadened fitburst model. + + Parameters + ---------- + burst_width : float + the temporal width of the burst component + + time_diff : float + the dispersion-delayed timeseries relative to the arrival time + + scattering_timescale : float + the scattering timescale at the frequency corresponding to 'time_diff' + + Returns + ------- + arg_erf : float + the argument of the error function + """ + + # compute and return the argument of the error function. + arg_erf = (time_diff - burst_width ** 2 / scattering_timescale) / burst_width / np.sqrt(2) + + return arg_erf + +def argument_exp(burst_width: float, time_diff: float, scattering_timescale: float) -> float: + """ + Computes the argument of the exponential term that is common to all derivatives + of the scatter-broadened fitburst model. + + Parameters + ---------- + + burst_width : float + the temporal width of the burst component + + time_diff : float + the dispersion-delayed timeseries relative to the arrival time + + scattering_timescale : float + the scattering timescale at the frequency corresponding to 'time_diff' + + Returns + ------- + arg_exp : float + the argument of the exponential term + """ + + # compute the argument to the error function. + arg_erf = argument_erf(burst_width, time_diff, scattering_timescale) + + # compute and return the argument of the exponential term. + arg_exp = 0.5 * (burst_width / scattering_timescale) ** 2 - time_diff / scattering_timescale - arg_erf ** 2 + + return arg_exp + +def deriv_time_dm(name: str, freq: float, ref_freq: float, dm: float, dm_index: float): + """ + Computes the derivative of the dispersed timeseries with respect to either the DM + of DM-index parameters. + + Parameters + ---------- + name : str + name of the fit parameter for which to compute the partial derivative + + freq : float + value of the electromagnetic frequency at which to evaluate dispersed timeseries + + ref_freq : float + value of the electromagnetic frequency at which to reference dispersion + + dm : float + value of dispersion measure + + dm_index : float + value of exponent for frequency dependence of dispersion relation + + Returns + ------- + deriv_partial : float + partial derivative of the dispersed timeseries with respect to the 'name' fit parameter + """ + + deriv_partial = 0. + dm_const = general["constants"]["dispersion"] + + if name == "dm": + diff_freq = freq ** dm_index - ref_freq ** dm_index + deriv_partial = -dm_const * diff_freq + + elif name == "dm_index": + diff_freq = np.log(freq) * freq ** dm_index - np.log(ref_freq) * ref_freq ** dm_index + deriv_partial = -dm_const * dm * diff_freq + + return deriv_partial + +def deriv_argument_erf(name: str, freq: float, time_diff: float, parameters: dict, + component: int = 0) -> float: + """ + Computes the derivative of the error-function argument. + + Parameters + ---------- + name : str + name of the fit parameter for which to compute the partial derivative + + burst_width : float + the temporal width of the burst component + + time_diff : float + the dispersion-delayed timeseries relative to the arrival time + + sc_time : float + the scattering timescale at the frequency corresponding to 'time_diff' + + sc_time_ref : float + the scattering timescale at the frequency corresponding to 'time_diff' + + Returns + ------- + arg_exp : float + the argument of the exponential term + """ + + # get parameters. + width = parameters["burst_width"][component] + dm = parameters["dm"][0] + dm_index = parameters["dm_index"][0] + ref_freq = parameters["ref_freq"][component] + sc_time = parameters["scattering_timescale"][0] + sc_index = parameters["scattering_index"][0] + sc_time_freq = sc_time * (freq / ref_freq) ** sc_index + + # compute derivative with respect to the appropriate parameter. + deriv_first = 0. + + if name == "arrival_time": + deriv_first = -1 / width / np.sqrt(2) + + elif name == "burst_width": + deriv_first = -(time_diff / width ** 2 + 1 / sc_time_freq) / np.sqrt(2) + + elif name == "dm" or name == "dm_index": + deriv_first = deriv_time_dm(name, freq, ref_freq, dm, dm_index) / width / np.sqrt(2) + + elif name == "scattering_timescale": + deriv_first = width / sc_time / sc_time_freq / np.sqrt(2) + + elif name == "scattering_index": + deriv_first = np.log(freq / ref_freq) * width / sc_time_freq / np.sqrt(2) + + return deriv_first + +def deriv_argument_exp(name: str, freq: float, time_diff: float, parameters: dict, + component: int = 0) -> float: + """ + Computes the argument of the exponential term that is common to all derivatives + of the scatter-broadened fitburst model. + + Parameters + ---------- + name : str + name of the fit parameter for which to compute the partial derivative + + burst_width : float + the temporal width of the burst component + + time_diff : float + the dispersion-delayed timeseries relative to the arrival time + + sc_time : float + the scattering timescale at the frequency corresponding to 'time_diff' + + sc_time_ref : float + the scattering timescale at the frequency corresponding to 'time_diff' + + Returns + ------- + arg_exp : float + the argument of the exponential term + """ + + # get parameters. + width = parameters["burst_width"][component] + dm = parameters["dm"][0] + dm_index = parameters["dm_index"][0] + ref_freq = parameters["ref_freq"][component] + sc_time = parameters["scattering_timescale"][0] + sc_index = parameters["scattering_index"][0] + sc_time_freq = sc_time * (freq / ref_freq) ** sc_index + arg_erf = argument_erf(width, time_diff, sc_time_freq) + deriv_arg_erf = deriv_argument_erf(name, freq, time_diff, parameters, component) + + # compute derivative with respect to the appropriate parameter. + deriv_first = 0. + + if name == "arrival_time": + deriv_first = 1 / sc_time_freq + + elif name == "burst_width": + deriv_first = width / sc_time_freq ** 2 + + elif name == "dm" or name == "dm_index": + deriv_first = -deriv_time_dm(name, freq, ref_freq, dm, dm_index) / sc_time_freq + + elif name == "scattering_timescale": + deriv_first = -(width / sc_time_freq) ** 2 / sc_time + time_diff / sc_time_freq / sc_time + + elif name == "scattering_index": + log_freq = np.log(freq / ref_freq) + deriv_first = log_freq * (-(width / sc_time_freq) ** 2 + time_diff / sc_time_freq) + + # now apply subtraction that is common to all derivatives. + deriv_first -= 2 * arg_erf * deriv_arg_erf + + return deriv_first + +def deriv2_argument_erf(name1: str, name2: str, freq: float, time_diff: float, parameters: dict, + component: int = 0) -> float: + """ + Computes the mixed-partial derivative of the error-function argument with + respect to one or two fitburst parameters that define a scatter-broadened model. + + Parameters + ---------- + name1 : str + name of the fit parameter for which to compute the partial derivative + + name2 : str + name of the fit parameter for which to compute the partial derivative + + burst_width : float + the temporal width of the burst component + + time_diff : float + the dispersion-delayed timeseries relative to the arrival time + + sc_time : float + the scattering timescale at the frequency corresponding to 'time_diff' + + sc_time_ref : float + the scattering timescale at the frequency corresponding to 'time_diff' + + Returns + ------- + arg_exp : float + the argument of the exponential term + """ + + # get parameters. + width = parameters["burst_width"][component] + dm = parameters["dm"][0] + dm_index = parameters["dm_index"][0] + ref_freq = parameters["ref_freq"][component] + sc_time = parameters["scattering_timescale"][0] + sc_index = parameters["scattering_index"][0] + sc_time_freq = sc_time * (freq / ref_freq) ** sc_index + deriv_mixed = 0. + + if name1 == "arrival_time" and name2 == "burst_width": + deriv_mixed = 1 / width ** 2 / np.sqrt(2) + + elif name1 == "burst_width" and name2 == "burst_width": + deriv_mixed = np.sqrt(2) * time_diff / width ** 3 + + elif (name1 == "dm" and name2 == "burst_width") or (name1 == "dm_index" and name2 == "burst_width"): + deriv_mixed = -deriv_time_dm(name1, freq, ref_freq, dm, dm_index) / width ** 2 / np.sqrt(2) + + elif name1 == "scattering_index" and name2 == "burst_width": + log_freq = np.log(freq / ref_freq) + deriv_mixed = log_freq / np.sqrt(2) / sc_time_freq + + elif name1 == "scattering_timescale" and name2 == "burst_width": + deriv_mixed = 1 / sc_time / sc_time_freq / np.sqrt(2) + + return deriv_mixed + def deriv_model_wrt_amplitude(parameters: dict, model: object, component: int = 0) -> float: """ Computes the derivative of the model with respect to the amplitude parameter. @@ -133,14 +413,16 @@ def deriv_model_wrt_burst_width(parameters: dict, model: float, component: int = # define argument of error and scattering timescales over frequency. log_freq = np.log(freq_ratio[freq]) spectrum = 10 ** amplitude * freq_ratio[freq] ** (spectral_index + spectral_running * log_freq) + spectrum *= freq_ratio[freq] ** (-scattering_index) current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - erf_arg = (current_timediff - burst_width ** 2 / scat_times_freq[freq]) / burst_width / np.sqrt(2) - exp_arg = (burst_width / scat_times_freq[freq]) ** 2 / 2 - current_timediff / scat_times_freq[freq] - erf_arg ** 2 + arg_exp = argument_exp(burst_width, current_timediff, scat_times_freq[freq]) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff, parameters, component + ) # now compute derivative contribution from current component. - term1 = burst_width * model.spectrum_per_component[freq, :, component] / scat_times_freq[freq] ** 2 - term2 = -np.sqrt(2 / np.pi) * spectrum * current_timediff * np.exp(exp_arg) / \ - burst_width ** 2 / scat_times_freq[freq] + term1 = (burst_width / scat_times_freq[freq] ** 2) * model.spectrum_per_component[freq, :, component] + term2 = spectrum * np.exp(arg_exp) * deriv_arg_erf * 2 / np.sqrt(np.pi) deriv_mod[freq, :] += term1 + term2 @@ -192,13 +474,16 @@ def deriv_model_wrt_arrival_time(parameters: dict, model: float, component: int # define argument of error and scattering timescales over frequency.a log_freq = np.log(freq_ratio[freq]) spectrum = 10 ** amplitude * freq_ratio[freq] ** (spectral_index + spectral_running * log_freq) + spectrum *= freq_ratio[freq] ** (-scattering_index) current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - erf_arg = (current_timediff - burst_width ** 2 / scat_times_freq[freq]) / burst_width / np.sqrt(2) - exp_arg = (burst_width / scat_times_freq[freq]) ** 2 / 2 - current_timediff / scat_times_freq[freq] - erf_arg ** 2 + arg_exp = argument_exp(burst_width, current_timediff, scat_times_freq[freq]) + deriv_arg_erf = deriv_argument_erf( + "arrival_time", model.freqs[freq], current_timediff, parameters, component + ) # now compute derivative contribution from current component. term1 = model.spectrum_per_component[freq, :, component] / scat_times_freq[freq] - term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width / scat_times_freq[freq] + term2 = 2 * spectrum * np.exp(arg_exp) * deriv_arg_erf / np.sqrt(np.pi) deriv_mod[freq, :] += term1 + term2 @@ -227,7 +512,7 @@ def deriv_model_wrt_dm(parameters: dict, model: float, component: int = 0, add_a # get dimensions and define empty model-derivative matrix. num_freq, num_time, num_component = model.timediff_per_component.shape deriv_mod_int = np.zeros((num_freq, num_time, num_component), dtype=float) - dm_const = 4149.377593360996 + dm = parameters["dm"][0] dm_index = parameters["dm_index"][0] scattering_index = parameters["scattering_index"][0] scattering_timescale = parameters["scattering_timescale"][0] @@ -237,35 +522,37 @@ def deriv_model_wrt_dm(parameters: dict, model: float, component: int = 0, add_a amplitude = parameters["amplitude"][current_component] burst_width = parameters["burst_width"][current_component] ref_freq = parameters["ref_freq"][current_component] + freq_ratio = model.freqs / ref_freq + sc_times_freq = scattering_timescale * freq_ratio ** scattering_index spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] timediff = model.timediff_per_component[:, :, current_component] - freq_ratio = model.freqs / ref_freq - freq_diff = model.freqs ** dm_index - ref_freq ** dm_index - scat_times_freq = scattering_timescale * freq_ratio ** scattering_index + deriv_timediff_wrt_dm = deriv_time_dm("dm", model.freqs, ref_freq, dm, dm_index) - for freq_idx in range(num_freq): - current_timediff = timediff[freq_idx, :] + for freq in range(num_freq): + current_timediff = timediff[freq, :] + current_model = model.spectrum_per_component[freq, :, current_component] - if scat_times_freq[freq_idx] < np.fabs(0.15 * burst_width): - deriv_timediff_wrt_dm = -dm_const * freq_diff[freq_idx] - deriv_mod_int[freq_idx, :, current_component] += (-current_timediff * - model.spectrum_per_component[freq_idx, :, current_component] / burst_width ** 2 * deriv_timediff_wrt_dm) + if sc_times_freq[freq] < np.fabs(0.15 * burst_width): + deriv_mod_int[freq, :, current_component] = (-current_timediff * + model.spectrum_per_component[freq, :, current_component] / burst_width ** 2 * deriv_timediff_wrt_dm[freq]) else: # define argument of error and scattering timescales over frequency. - log_freq = np.log(freq_ratio[freq_idx]) - spectrum = 10 ** amplitude * freq_ratio[freq_idx] ** (spectral_index + spectral_running * log_freq) + log_freq = np.log(freq_ratio[freq]) + spectrum = 10 ** amplitude * freq_ratio[freq] ** (spectral_index + spectral_running * log_freq) + spectrum *= freq_ratio[freq] ** (-scattering_index) current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - erf_arg = (current_timediff - burst_width ** 2 / scat_times_freq[freq_idx]) / burst_width / np.sqrt(2) - exp_arg = (burst_width / scat_times_freq[freq_idx]) ** 2 / 2 - current_timediff / scat_times_freq[freq_idx] - \ - erf_arg ** 2 + arg_exp = argument_exp(burst_width, current_timediff, sc_times_freq[freq]) + deriv_arg_erf = deriv_argument_erf( + "dm", model.freqs[freq], current_timediff, parameters, current_component + ) # now compute derivative contribution from current component. - term1 = model.spectrum_per_component[freq_idx, :, current_component] / scat_times_freq[freq_idx] - term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width / scat_times_freq[freq_idx] - deriv_mod_int[freq_idx, :, current_component] += dm_const * freq_diff[freq_idx] * (term1 + term2) + term1 = -deriv_timediff_wrt_dm[freq] * current_model / sc_times_freq[freq] + term2 = spectrum * np.exp(arg_exp) * deriv_arg_erf * 2 / np.sqrt(np.pi) + deriv_mod_int[freq, :, current_component] = (term1 + term2) # now determine if all components should be summed or not. deriv_mod = deriv_mod_int[:, :, component] @@ -296,36 +583,49 @@ def deriv_model_wrt_dm_index(parameters: dict, model: float, component: int = 0, """ # get dimensions and define empty model-derivative matrix. + num_freq, num_time, num_component = model.timediff_per_component.shape deriv_mod_int = np.zeros(model.spectrum_per_component.shape) dm = parameters["dm"][0] - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] sc_index = parameters["scattering_index"][0] sc_time = parameters["scattering_timescale"][0] - for current_component in range(model.spectrum_per_component.shape[2]): + for current_component in range(num_component): amplitude = parameters["amplitude"][current_component] burst_width = parameters["burst_width"][current_component] - current_timediff = model.timediff_per_component[:, : , current_component] + timediff = model.timediff_per_component[:, : , current_component] ref_freq = parameters["ref_freq"][current_component] + freq_ratio = model.freqs / ref_freq spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] + sc_times_freq = sc_time * freq_ratio ** sc_index + deriv_timediff_wrt_dm = deriv_time_dm("dm_index", model.freqs, ref_freq, dm, dm_index) + for freq in range (num_freq): + current_timediff = timediff[freq, :] - for freq_idx in range (model.spectrum_per_component.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq - sc_time_freq = sc_time * freq_ratio ** sc_index - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - product = np.log(model.freqs[freq_idx]) * model.freqs[freq_idx] ** dm_index -\ - np.log(ref_freq) * ref_freq ** dm_index + if sc_times_freq[freq] < np.fabs(0.15 * burst_width): + deriv_mod_int[freq, :, current_component] += (-current_timediff * + model.spectrum_per_component[freq, :, current_component] / burst_width ** 2 * deriv_timediff_wrt_dm[freq]) + + else: - term1 = dm_const * dm * product * \ - model.spectrum_per_component[freq_idx, :, current_component] / sc_time_freq - term2 = -np.sqrt(2 / np.pi) * product * spectrum * np.exp(exp_arg) / burst_width / sc_time_freq + # define argument of error and scattering timescales over frequency. + log_freq = np.log(freq_ratio[freq]) + spectrum = 10 ** amplitude * freq_ratio[freq] ** (spectral_index + spectral_running * log_freq) + spectrum *= freq_ratio[freq] ** (-sc_index) + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width + arg_exp = argument_exp(burst_width, current_timediff, sc_times_freq[freq]) + deriv_arg_erf = deriv_argument_erf( + "dm_index", model.freqs[freq], current_timediff, parameters, current_component + ) + + # now calculate terms in partial derivative. + term1 = -deriv_timediff_wrt_dm[freq] * model.spectrum_per_component[freq, :, current_component] + term1 /= sc_times_freq[freq] + term2 = spectrum * np.exp(arg_exp) * deriv_arg_erf * 2 / np.sqrt(np.pi) - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + deriv_mod_int[freq, :, current_component] = (term1 + term2) # now determine if all components should be summed or not. deriv_mod = deriv_mod_int[:, :, component] @@ -368,29 +668,30 @@ def deriv_model_wrt_scattering_timescale(parameters: dict, model: float, compone amplitude = parameters["amplitude"][current_component] burst_width = parameters["burst_width"][current_component] ref_freq = parameters["ref_freq"][current_component] + freq_ratio = model.freqs / ref_freq + log_freq = np.log(freq_ratio) scattering_index = parameters["scattering_index"][0] scattering_timescale = parameters["scattering_timescale"][0] spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - timediff = model.timediff_per_component[:, :, current_component] + current_timediff = model.timediff_per_component[:, :, current_component] # define argument of error and scattering timescales over frequency. - freq_ratio = model.freqs / ref_freq - log_freq = np.log(freq_ratio) spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * log_freq) - timediff[timediff < -5 * burst_width] = -5 * burst_width + spectrum *= freq_ratio ** (-scattering_index) + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width scat_times_freq = scattering_timescale * freq_ratio ** scattering_index - erf_arg = (timediff - burst_width ** 2 / scat_times_freq[:, None]) / burst_width / np.sqrt(2) - product = np.exp((burst_width / scat_times_freq[:, None]) ** 2 / 2 - timediff / scat_times_freq[:, None] - - erf_arg ** 2) + arg_exp = argument_exp(burst_width, current_timediff, scat_times_freq[:, None]) + deriv_arg_erf = deriv_argument_erf( + "scattering_timescale", model.freqs[:, None], current_timediff, parameters, current_component + ) # now compute derivative contribution from current component. - term1 = model.spectrum_per_component[:, :, current_component] - term2 = (burst_width / scat_times_freq[:, None]) ** 2 * model.spectrum_per_component[:, :, current_component] - term3 = -timediff / scat_times_freq[:, None] * model.spectrum_per_component[:, :, current_component] - term4 = -np.sqrt(2 / np.pi) * product * burst_width / scat_times_freq[:, None] ** 2 + term1 = (-(burst_width / scat_times_freq[:, None]) ** 2 + current_timediff / scat_times_freq[:, None]) * \ + model.spectrum_per_component[:, :, current_component] / scattering_timescale + term2 = spectrum[:, None] * np.exp(arg_exp) * deriv_arg_erf * 2 / np.sqrt(np.pi) - deriv_mod_int[:, :, current_component] = -(term1 + term2 + term3 + term4 * spectrum[:, None]) / scattering_timescale + deriv_mod_int[:, :, current_component] = term1 + term2 # now determine if all components should be summed or not. deriv_mod = deriv_mod_int[:, :, component] @@ -433,29 +734,30 @@ def deriv_model_wrt_scattering_index(parameters: dict, model: float, component: amplitude = parameters["amplitude"][current_component] burst_width = parameters["burst_width"][current_component] ref_freq = parameters["ref_freq"][current_component] + freq_ratio = model.freqs / ref_freq + log_freq = np.log(freq_ratio) scattering_index = parameters["scattering_index"][0] scattering_timescale = parameters["scattering_timescale"][0] spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - timediff = model.timediff_per_component[:, :, current_component] + current_timediff = model.timediff_per_component[:, :, current_component] + scat_times_freq = scattering_timescale * freq_ratio ** scattering_index + spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * log_freq) + spectrum *= freq_ratio ** (-scattering_index) # define argument of error and scattering timescales over frequency. - freq_ratio = model.freqs / ref_freq - log_freq = np.log(freq_ratio) - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * log_freq) - timediff[timediff < -5 * burst_width] = -5 * burst_width - scat_times_freq = scattering_timescale * freq_ratio ** scattering_index - erf_arg = (timediff - burst_width ** 2 / scat_times_freq[:, None]) / burst_width / np.sqrt(2) - product = np.exp((burst_width / scat_times_freq[:, None]) ** 2 / 2 - timediff / scat_times_freq[:, None] - - erf_arg ** 2) + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width + arg_exp = argument_exp(burst_width, current_timediff, scat_times_freq[:, None]) + deriv_arg_erf = deriv_argument_erf( + "scattering_index", model.freqs, current_timediff, parameters, current_component + ) # now compute derivative contribution from current component. - term1 = model.spectrum_per_component[:, :, current_component] - term2 = (burst_width / scat_times_freq[:, None]) ** 2 * model.spectrum_per_component[:, :, current_component] - term3 = -timediff / scat_times_freq[:, None] * model.spectrum_per_component[:, :, current_component] - term4 = -np.sqrt(2 / np.pi) * product * burst_width / scat_times_freq[:, None] ** 2 + term1 = -log_freq * (1 + (burst_width / scat_times_freq[:, None]) ** 2 - current_timediff / scat_times_freq[:, None]) + term1 *= model.spectrum_per_component[:, :, current_component] + term2 = spectrum * np.exp(arg_exp) * deriv_arg_erf[:, None] * 2 / np.sqrt(np.pi) - deriv_mod_int[:, :, current_component] = -(term1 + term2 + term3 + term4 * spectrum[:, None]) * log_freq[:, None] + deriv_mod_int[:, :, current_component] = (term1 + term2) # now determine if all components should be summed or not. deriv_mod = deriv_mod_int[:, :, component] @@ -1134,40 +1436,41 @@ def deriv2_model_wrt_burst_width_burst_width(parameters: dict, model: float, com sc_time = parameters["scattering_timescale"][0] # global parameter. spectral_index = parameters["spectral_index"][component] spectral_running = parameters["spectral_running"][component] + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index # if scattering is not resolvable, then assume Gaussian temporal profile. if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod[freq_idx, :] = -3 * current_timediff[freq_idx, :] ** 2 * \ - current_model[freq_idx, :] / burst_width ** 4 - deriv_mod[freq_idx, :] += (current_timediff[freq_idx, :] ** 4 * \ - current_model[freq_idx, :] / burst_width ** 6) + deriv_mod[freq, :] = -3 * current_timediff[freq, :] ** 2 * \ + current_model[freq, :] / burst_width ** 4 + deriv_mod[freq, :] += (current_timediff[freq, :] ** 4 * \ + current_model[freq, :] / burst_width ** 6) else: # adjust time-difference values to make them friendly for error function. - current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = -(current_timediff[freq_idx, :] / burst_width ** 2 + 1 / sc_time_freq) / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = burst_width / sc_time_freq ** 2 - 2 * erf_arg * erf_arg_deriv + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "burst_width", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contribute to mixed derivative. - term1 = current_model[freq_idx, :] / sc_time_freq ** 2 - term2 = burst_width * deriv_first[freq_idx, :] / sc_time_freq ** 2 - term3 = 2 * np.sqrt(2) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) /\ - sc_time_freq / burst_width ** 3 / np.sqrt(np.pi) - term4 = -np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) *\ - exp_arg_deriv / sc_time_freq / burst_width ** 2 - - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term4 + term1 = current_model[freq, :] / sc_time_freq ** 2 + term2 = burst_width * deriv_first[freq, :] / sc_time_freq ** 2 + term3 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term4 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod[freq, :] = term1 + term2 + term3 + term4 return deriv_mod @@ -1206,33 +1509,36 @@ def deriv2_model_wrt_burst_width_arrival_time(parameters: dict, model: float, co spectral_running = parameters["spectral_running"][component] # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index # if scattering is not resolvable, then assume Gaussian temporal profile. if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod[freq_idx, :] = -2 * current_timediff[freq_idx, :] * current_model[freq_idx, :] / burst_width ** 3 - deriv_mod[freq_idx, :] += (current_timediff[freq_idx, :] ** 3 * current_model[freq_idx, :] / burst_width ** 5) + deriv_mod[freq, :] = -2 * current_timediff[freq, :] * current_model[freq, :] / burst_width ** 3 + deriv_mod[freq, :] += (current_timediff[freq, :] ** 3 * current_model[freq, :] / burst_width ** 5) else: # adjust time-difference values to make them friendly for error function. current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = -1 / burst_width / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = 1 / sc_time_freq - 2 * erf_arg * erf_arg_deriv + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "arrival_time", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "arrival_time", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. - term1 = burst_width * deriv_first[freq_idx, :] / sc_time_freq ** 2 - term2 = np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width ** 2 / sc_time_freq - term3 = -np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) *\ - exp_arg_deriv / burst_width ** 2 / sc_time_freq - - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term1 = burst_width * deriv_first[freq, :] / sc_time_freq ** 2 + term2 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term3 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod[freq, :] = term1 + term2 + term3 return deriv_mod @@ -1274,28 +1580,28 @@ def deriv2_model_wrt_burst_width_scattering_timescale(parameters: dict, model: f current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_ratio_sc = freq_ratio ** sc_index - sc_time_freq = sc_time * freq_ratio_sc + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq + sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = burst_width * freq_ratio_sc / sc_time_freq ** 2 / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = -burst_width ** 2 * freq_ratio_sc / sc_time_freq ** 3 + \ - current_timediff[freq_idx, :] / sc_time_freq ** 2 * freq_ratio_sc - \ - 2 * erf_arg * erf_arg_deriv + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "scattering_timescale", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. - term1 = -2 * burst_width * current_model[freq_idx, :] * freq_ratio_sc / sc_time_freq ** 3 - term2 = burst_width * deriv_first[freq_idx, :] / sc_time_freq ** 2 - term3 = np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * freq_ratio_sc * spectrum * np.exp(exp_arg) /\ - burst_width ** 2 / sc_time_freq ** 2 - term4 = -np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) * exp_arg_deriv /\ - burst_width ** 2 / sc_time_freq - - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term4 + term1 = -2 * burst_width * current_model[freq, :] * sc_time / sc_time_freq ** 2 + term2 = burst_width * deriv_first[freq, :] / sc_time_freq ** 2 + term3 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term4 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod[freq, :] = term1 + term2 + term3 + term4 return deriv_mod @@ -1337,28 +1643,31 @@ def deriv2_model_wrt_burst_width_scattering_index(parameters: dict, model: float current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_ratio_sc = freq_ratio ** sc_index - sc_time_freq = sc_time * freq_ratio_sc + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq + log_freq = np.log_freq_ratio + sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = burst_width * np.log(freq_ratio) / sc_time_freq / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = -(burst_width / sc_time_freq) ** 2 * np.log(freq_ratio) + \ - current_timediff[freq_idx, :] / sc_time_freq * np.log(freq_ratio) - 2 * erf_arg * erf_arg_deriv + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "scattering_index", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "scattering_index", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. - term1 = -2 * burst_width * current_model[freq_idx, :] * np.log(freq_ratio) / sc_time_freq ** 2 - term2 = burst_width * deriv_first[freq_idx, :] / sc_time_freq ** 2 - term3 = np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * np.log(freq_ratio) * spectrum * np.exp(exp_arg) /\ - burst_width ** 2 / sc_time_freq - term4 = -np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) * exp_arg_deriv /\ - burst_width ** 2 / sc_time_freq + term1 = -2 * log_freq * burst_width * current_model[freq, :] / sc_time_freq ** 2 + term2 = burst_width * deriv_first[freq, :] / sc_time_freq ** 2 + term3 = -spectrum * np.exp(arg_exp) * deriv_arg_erf * log_freq * 2 / np.sqrt(np.pi) + term4 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term5 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term4 + deriv_mod[freq, :] = term1 + term2 + term3 + term4 + term5 return deriv_mod @@ -1387,9 +1696,10 @@ def deriv2_model_wrt_burst_width_dm(parameters: dict, model: float, component: i burst_width = parameters["burst_width"][component] current_model = model.spectrum_per_component[:, :, component] current_timediff = model.timediff_per_component[:, :, component] + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width deriv_first = deriv_model_wrt_dm(parameters, model, component, add_all = False) deriv_mod = np.zeros(current_model.shape) - dm_const = 4149.377593360996 + dm = parameters["dm"][0] # global parameter. dm_index = parameters["dm_index"][0] # global parameter. ref_freq = parameters["ref_freq"][component] sc_index = parameters["scattering_index"][0] # global parameter. @@ -1398,35 +1708,37 @@ def deriv2_model_wrt_burst_width_dm(parameters: dict, model: float, component: i spectral_running = parameters["spectral_running"][component] # global parameter. # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq - product = dm_const * (model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index) + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index # if scattering is not resolvable, then assume Gaussian temporal profile. if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod[freq_idx, :] += -2 * current_timediff[freq_idx, :] * product * \ - current_model[freq_idx, :] / burst_width ** 3 - deriv_mod[freq_idx, :] += (current_timediff[freq_idx, :] ** 2 * deriv_first[freq_idx, :] / burst_width ** 3) + deriv_time_dm_freq = deriv_time_dm("dm", model.freqs[freq], ref_freq, dm, dm_index) + deriv_mod[freq, :] += -2 * current_timediff[freq, :] * deriv_time_dm_freq * \ + current_model[freq, :] / burst_width ** 3 + deriv_mod[freq, :] += (current_timediff[freq, :] ** 2 * deriv_first[freq, :] / burst_width ** 3) else: # adjust time-difference values to make them friendly for error function. - current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = -product / burst_width / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = product / sc_time_freq - 2 * erf_arg * erf_arg_deriv + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "dm", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "dm", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. - term1 = burst_width * deriv_first[freq_idx, :] / sc_time_freq ** 2 - term2 = np.sqrt(2 / np.pi) * spectrum * product * np.exp(exp_arg) / burst_width ** 2 / sc_time_freq - term3 = -np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) *\ - exp_arg_deriv / burst_width ** 2 / sc_time_freq - - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term1 = burst_width * current_model[freq, :] / sc_time_freq ** 2 + term2 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term3 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod[freq, :] = term1 + term2 + term3 return deriv_mod @@ -1455,10 +1767,10 @@ def deriv2_model_wrt_burst_width_dm_index(parameters: dict, model: float, compon burst_width = parameters["burst_width"][component] current_model = model.spectrum_per_component[:, :, component] current_timediff = model.timediff_per_component[:, :, component] + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width deriv_first = deriv_model_wrt_dm_index(parameters, model, component, add_all = False) deriv_mod = np.zeros(current_model.shape) dm = parameters["dm"][0] # global parameter. - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] # global parameter. ref_freq = parameters["ref_freq"][component] dm_index = parameters["dm_index"][0] # global parameter. @@ -1467,36 +1779,38 @@ def deriv2_model_wrt_burst_width_dm_index(parameters: dict, model: float, compon spectral_index = parameters["spectral_index"][component] # global parameter. spectral_running = parameters["spectral_running"][component] # global parameter. - # adjust time-difference values to make them friendly for error function. - current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq - product = dm_const * dm * (np.log(model.freqs[freq_idx]) * model.freqs[freq_idx] ** dm_index - \ - np.log(ref_freq) * ref_freq ** dm_index) + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index # if scattering is not resolvable, then assume Gaussian temporal profile. if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod[freq_idx, :] += -2 * current_timediff[freq_idx, :] * product * current_model[freq_idx, :] / burst_width ** 3 - deriv_mod[freq_idx, :] += (current_timediff[freq_idx, :] ** 2 * deriv_first[freq_idx, :] / burst_width ** 3) + deriv_time_dm_freq = deriv_time_dm("dm_index", model.freqs[freq], ref_freq, dm, dm_index) + deriv_mod[freq, :] += -2 * current_timediff[freq, :] * deriv_time_dm_freq * current_model[freq, :] / burst_width ** 3 + deriv_mod[freq, :] += (current_timediff[freq, :] ** 2 * deriv_first[freq, :] / burst_width ** 3) else: + # adjust time-difference values to make them friendly for error function. spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = -product / burst_width / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = product / sc_time_freq - 2 * erf_arg * erf_arg_deriv + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "dm_index", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "dm_index", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. - term1 = burst_width * deriv_first[freq_idx, :] / sc_time_freq ** 2 - term2 = np.sqrt(2 / np.pi) * spectrum * product * np.exp(exp_arg) / burst_width ** 2 / sc_time_freq - term3 = -np.sqrt(2 / np.pi) * current_timediff[freq_idx, :] * spectrum * np.exp(exp_arg) *\ - exp_arg_deriv / burst_width ** 2 / sc_time_freq + term1 = burst_width * deriv_first[freq, :] / sc_time_freq ** 2 + term2 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term3 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) - deriv_mod[freq_idx, :] = term1 + term2 + term3 + deriv_mod[freq, :] = term1 + term2 + term3 return deriv_mod @@ -1534,24 +1848,25 @@ def deriv2_model_wrt_arrival_time_arrival_time(parameters: dict, model: float, c deriv_mod = np.zeros((num_freq, num_time), dtype=float) - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod[freq_idx, :] = current_timediff[freq_idx, :] ** 2 * current_model[freq_idx, :] / burst_width ** 4 - deriv_mod[freq_idx, :] -= (current_model[freq_idx, :] / burst_width ** 2) + deriv_mod[freq, :] = current_timediff[freq, :] ** 2 * current_model[freq, :] / burst_width ** 4 + deriv_mod[freq, :] -= (current_model[freq, :] / burst_width ** 2) else: spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + spectrum *= freq_ratio ** (-sc_index) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = -1 / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 exp_arg_deriv = (1 / sc_time_freq - 2 * erf_arg * erf_arg_deriv) - term1 = deriv_first[freq_idx, :] / sc_time_freq - term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width / sc_time_freq - deriv_mod[freq_idx, :] = term1 + term2 + term1 = deriv_first[freq, :] / sc_time_freq + term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + deriv_mod[freq, :] = term1 + term2 return deriv_mod @@ -1591,31 +1906,32 @@ def deriv2_model_wrt_arrival_time_dm(parameters: dict, model: float, component: spectral_running = parameters["spectral_running"][component] # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod[freq_idx, :] = current_timediff[freq_idx, :] * deriv_first[freq_idx, :] / burst_width ** 2 - deriv_mod[freq_idx, :] -= (dm_const * (model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index) * - current_model[freq_idx, :] / burst_width ** 2) + deriv_mod[freq, :] = current_timediff[freq, :] * deriv_first[freq, :] / burst_width ** 2 + deriv_mod[freq, :] -= (dm_const * (model.freqs[freq] ** dm_index - ref_freq ** dm_index) * + current_model[freq, :] / burst_width ** 2) else: # adjust time-difference values to make them friendly for error function. current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = -dm_const * (model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index) / burst_width / np.sqrt(2) + spectrum *= freq_ratio ** (-sc_index) + erf_arg = (current_timediff[freq, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) + erf_arg_deriv = -dm_const * (model.freqs[freq] ** dm_index - ref_freq ** dm_index) / burst_width / np.sqrt(2) exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = dm_const * (model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index) / sc_time_freq - \ + current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 + exp_arg_deriv = dm_const * (model.freqs[freq] ** dm_index - ref_freq ** dm_index) / sc_time_freq - \ 2 * erf_arg * erf_arg_deriv # now define terms that contriubte to mixed-partial derivative. - term1 = deriv_first[freq_idx, :] / sc_time_freq - term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width / sc_time_freq - deriv_mod[freq_idx, :] = term1 + term2 + term1 = deriv_first[freq, :] / sc_time_freq + term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + deriv_mod[freq, :] = term1 + term2 return deriv_mod @@ -1695,24 +2011,25 @@ def deriv2_model_wrt_arrival_time_scattering_timescale(parameters: dict, model: current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) + spectrum *= freq_ratio ** (-sc_index) + erf_arg = (current_timediff[freq, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) erf_arg_deriv = burst_width * freq_ratio ** sc_index / sc_time_freq ** 2 / np.sqrt(2) exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 exp_arg_deriv = -burst_width ** 2 * freq_ratio ** sc_index / sc_time_freq ** 3 + \ - current_timediff[freq_idx, :] * freq_ratio ** sc_index / sc_time_freq ** 2 - \ + current_timediff[freq, :] * freq_ratio ** sc_index / sc_time_freq ** 2 - \ 2 * erf_arg * erf_arg_deriv # now define terms that contriubte to mixed-partial derivative. - term1 = -model.spectrum_per_component[freq_idx, :, component] * freq_ratio ** sc_index / sc_time_freq ** 2 - term2 = deriv_first[freq_idx, :] / sc_time_freq - term3 = np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * freq_ratio ** sc_index / burst_width / sc_time_freq ** 2 - term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width / sc_time_freq - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term4 + term1 = -model.spectrum_per_component[freq, :, component] / sc_time / sc_time_freq + term2 = deriv_first[freq, :] / sc_time_freq + term3 = 0. # np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * freq_ratio ** sc_index / burst_width / sc_time_freq ** 2 + term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + deriv_mod[freq, :] = term1 + term2 + term3 + term4 return deriv_mod @@ -1757,24 +2074,25 @@ def deriv2_model_wrt_arrival_time_scattering_index(parameters: dict, model: floa current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) + spectrum *= freq_ratio ** (-sc_index) + erf_arg = (current_timediff[freq, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) erf_arg_deriv = burst_width * np.log(freq_ratio) / sc_time_freq / np.sqrt(2) exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 exp_arg_deriv = -burst_width ** 2 * np.log(freq_ratio) / sc_time_freq ** 2 + \ - current_timediff[freq_idx, :] * np.log(freq_ratio) / sc_time_freq - \ + current_timediff[freq, :] * np.log(freq_ratio) / sc_time_freq - \ 2 * erf_arg * erf_arg_deriv # now define terms that contriubte to mixed-partial derivative. - term1 = -model.spectrum_per_component[freq_idx, :, component] * np.log(freq_ratio) / sc_time_freq - term2 = deriv_first[freq_idx, :] / sc_time_freq - term3 = np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * np.log(freq_ratio) / burst_width / sc_time_freq - term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width / sc_time_freq - deriv_mod[freq_idx, :] = term1 + term2 + term3 + term4 + term1 = -model.spectrum_per_component[freq, :, component] * np.log(freq_ratio) / sc_time_freq + term2 = deriv_first[freq, :] / sc_time_freq + term3 = np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * np.log(freq_ratio) / burst_width * freq_ratio ** (-sc_index) + term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + deriv_mod[freq, :] = term1 + term2 + term3 + term4 return deriv_mod @@ -1854,30 +2172,31 @@ def deriv2_model_wrt_dm_dm(parameters: dict, model: float, component: int = 0) - spectral_running = parameters["spectral_running"][current_component] # now loop over each frequency and compute mixed-derivative array per channel. - for freq_idx in range(current_model.shape[0]): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_diff = model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = model.freqs[freq] ** dm_index - ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod_int[freq_idx, :, current_component] = (dm_const * freq_diff) ** 2 * current_timediff[freq_idx, :] ** 2 * \ - current_model[freq_idx, :] / burst_width ** 4 - deriv_mod_int[freq_idx, :, current_component] -= (dm_const * freq_diff) ** 2 * current_model[freq_idx, :] / \ + deriv_mod_int[freq, :, current_component] = (dm_const * freq_diff) ** 2 * current_timediff[freq, :] ** 2 * \ + current_model[freq, :] / burst_width ** 4 + deriv_mod_int[freq, :, current_component] -= (dm_const * freq_diff) ** 2 * current_model[freq, :] / \ burst_width ** 2 else: # adjust time-difference values to make them friendly for error function. current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = -dm_const * freq_diff / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) + spectrum *= freq_ratio ** (-sc_index) exp_arg_deriv = dm_const * freq_diff / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = dm_const * freq_diff * deriv_first[freq_idx, :] / sc_time_freq - term2 = -dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width / \ - sc_time_freq * exp_arg_deriv - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term1 = dm_const * freq_diff * deriv_first[freq, :] / sc_time_freq + term2 = -dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width * \ + freq_ratio ** (-sc_index) * exp_arg_deriv + deriv_mod_int[freq, :, current_component] = term1 + term2 return np.sum(deriv_mod_int, axis = 2) @@ -1919,27 +2238,28 @@ def deriv2_model_wrt_scattering_timescale_scattering_timescale(parameters: dict, spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = (-1 / sc_time - burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq_idx, :] / \ + term0 = (-burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq, :] / \ sc_time / sc_time_freq) - term0_deriv = (1 / sc_time ** 2 + 3 * burst_width ** 2 / (sc_time * sc_time_freq) ** 2 - - 2 * current_timediff[freq_idx, :] / sc_time ** 2 / sc_time_freq) - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + term0_deriv = (3 * burst_width ** 2 / (sc_time * sc_time_freq) ** 2 - + 2 * current_timediff[freq, :] / sc_time ** 2 / sc_time_freq) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = burst_width / sc_time_freq / sc_time / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - exp_arg_deriv = -(burst_width / sc_time_freq )** 2 / sc_time + current_timediff[freq_idx, :] / sc_time_freq / sc_time - \ + spectrum *= freq_ratio ** (-sc_index) + exp_arg_deriv = -(burst_width / sc_time_freq )** 2 / sc_time + current_timediff[freq, :] / sc_time_freq / sc_time - \ 2 * erf_arg * erf_arg_deriv - term1 = term0_deriv * model.spectrum_per_component[freq_idx, :, current_component] - term2 = term0 * deriv_first[freq_idx, :] - term3 = -np.sqrt(8 / np.pi) * spectrum * burst_width * np.exp(exp_arg) / sc_time_freq / sc_time ** 2 - term4 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * exp_arg_deriv / sc_time_freq / sc_time + term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] + term2 = term0 * deriv_first[freq, :] + term3 = 0. # -np.sqrt(8 / np.pi) * spectrum * burst_width * np.exp(exp_arg) / sc_time_freq / sc_time ** 2 + term4 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * exp_arg_deriv / sc_time_freq - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + term4 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4 return np.sum(deriv_mod_int, axis = 2) @@ -1981,26 +2301,27 @@ def deriv2_model_wrt_scattering_timescale_scattering_index(parameters: dict, mod spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = (-1 / sc_time - burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq_idx, :] / \ + term0 = (-1 / sc_time - burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq, :] / \ sc_time / sc_time_freq) term0_deriv = 2 * (burst_width / sc_time_freq) ** 2 / sc_time * np.log(freq_ratio) - \ - current_timediff[freq_idx, :] * np.log(freq_ratio) / sc_time_freq / sc_time - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + current_timediff[freq, :] * np.log(freq_ratio) / sc_time_freq / sc_time + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = burst_width * np.log(freq_ratio) / sc_time_freq / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - exp_arg_deriv = -(burst_width / sc_time_freq) ** 2 * np.log(freq_ratio) + current_timediff[freq_idx, :] * \ + spectrum *= freq_ratio ** (-sc_index) + exp_arg_deriv = -(burst_width / sc_time_freq) ** 2 * np.log(freq_ratio) + current_timediff[freq, :] * \ np.log(freq_ratio) / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = term0_deriv * model.spectrum_per_component[freq_idx, :, current_component] - term2 = term0 * deriv_first[freq_idx, :] + term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] + term2 = term0 * deriv_first[freq, :] term3 = -np.sqrt(2 / np.pi) * spectrum * burst_width * np.log(freq_ratio) * np.exp(exp_arg) / sc_time_freq / sc_time term4 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * exp_arg_deriv / sc_time_freq / sc_time - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + term4 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4 return np.sum(deriv_mod_int, axis = 2) @@ -2044,23 +2365,24 @@ def deriv2_model_wrt_scattering_timescale_dm(parameters: dict, model: float, com spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_diff = model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = model.freqs[freq] ** dm_index - ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + spectrum *= freq_ratio ** (-sc_index) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = burst_width / sc_time / sc_time_freq / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = -1 / sc_time * (burst_width / sc_time_freq) ** 2 + current_timediff[freq_idx, :] / sc_time / sc_time_freq \ + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 + exp_arg_deriv = -1 / sc_time * (burst_width / sc_time_freq) ** 2 + current_timediff[freq, :] / sc_time / sc_time_freq \ -2 * erf_arg * erf_arg_deriv - term1 = dm_const * freq_diff * deriv_first[freq_idx, :] / sc_time_freq - term2 = -dm_const * freq_diff * current_model[freq_idx, :] / sc_time / sc_time_freq - term3 = dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width / sc_time_freq / sc_time - term4 = -dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width / sc_time_freq + term1 = dm_const * freq_diff * deriv_first[freq, :] / sc_time_freq + term2 = -dm_const * freq_diff * current_model[freq, :] / sc_time / sc_time_freq + term3 = dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width / freq_ratio ** (-sc_index) / sc_time + term4 = -dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + term4 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4 return np.sum(deriv_mod_int, axis = 2) @@ -2103,25 +2425,26 @@ def deriv2_model_wrt_scattering_timescale_dm_index(parameters: dict, model: floa spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_diff = np.log(model.freqs[freq_idx]) * model.freqs[freq_idx] ** dm_index - \ + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = np.log(model.freqs[freq]) * model.freqs[freq] ** dm_index - \ np.log(ref_freq) * ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = (-1 / sc_time - burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq_idx, :] / \ + term0 = (-1 / sc_time - burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq, :] / \ sc_time / sc_time_freq) term0_deriv = -dm_const * dm * freq_diff / sc_time_freq / sc_time - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = -dm_const * dm * freq_diff / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) + spectrum *= freq_ratio ** (-sc_index) exp_arg_deriv = dm_const * dm * freq_diff / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = term0_deriv * model.spectrum_per_component[freq_idx, :, current_component] - term2 = term0 * deriv_first[freq_idx, :] + term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] + term2 = term0 * deriv_first[freq, :] term3 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * exp_arg_deriv / sc_time_freq / sc_time - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 return np.sum(deriv_mod_int, axis = 2) @@ -2163,27 +2486,28 @@ def deriv2_model_wrt_scattering_index_scattering_index(parameters: dict, model: spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_diff = np.log(model.freqs[freq_idx]) * model.freqs[freq_idx] ** dm_index - \ + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = np.log(model.freqs[freq]) * model.freqs[freq] ** dm_index - \ np.log(ref_freq) * ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = -(1 + burst_width ** 2 / sc_time_freq ** 2 - current_timediff[freq_idx, :] / sc_time_freq) * np.log(freq_ratio) + term0 = -(1 + burst_width ** 2 / sc_time_freq ** 2 - current_timediff[freq, :] / sc_time_freq) * np.log(freq_ratio) term0_deriv = -np.log(freq_ratio) * (-2 * np.log(freq_ratio) * (burst_width / sc_time_freq) ** 2 + \ - current_timediff[freq_idx, :] * np.log(freq_ratio) / sc_time_freq) - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + current_timediff[freq, :] * np.log(freq_ratio) / sc_time_freq) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = burst_width * np.log(freq_ratio) / sc_time_freq / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) + spectrum *= freq_ratio ** (-sc_index) exp_arg_deriv = -np.log(freq_ratio) * (burst_width / sc_time_freq) ** 2 + \ - current_timediff[freq_idx, :] * np.log(freq_ratio) / sc_time_freq - 2 * erf_arg * erf_arg_deriv + current_timediff[freq, :] * np.log(freq_ratio) / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = term0_deriv * model.spectrum_per_component[freq_idx, :, current_component] - term2 = term0 * deriv_first[freq_idx, :] + term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] + term2 = term0 * deriv_first[freq, :] term3 = -np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * np.log(freq_ratio) ** 2 / sc_time_freq term4 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * np.log(freq_ratio) * exp_arg_deriv / sc_time_freq - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + term4 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4 return np.sum(deriv_mod_int, axis = 2) @@ -2226,22 +2550,23 @@ def deriv2_model_wrt_scattering_index_dm(parameters: dict, model: float, compone spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_diff = model.freqs[freq_idx] ** dm_index - ref_freq ** dm_index + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = model.freqs[freq] ** dm_index - ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = -(1 + burst_width ** 2 / sc_time_freq ** 2 - current_timediff[freq_idx, :] / sc_time_freq) * np.log(freq_ratio) + term0 = -(1 + burst_width ** 2 / sc_time_freq ** 2 - current_timediff[freq, :] / sc_time_freq) * np.log(freq_ratio) term0_deriv = -np.log(freq_ratio) * dm_const * freq_diff / sc_time_freq - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = -dm_const * freq_diff / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) + spectrum *= freq_ratio ** (-sc_index) exp_arg_deriv = dm_const * freq_diff / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = term0_deriv * model.spectrum_per_component[freq_idx, :, current_component] - term2 = term0 * deriv_first[freq_idx, :] + term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] + term2 = term0 * deriv_first[freq, :] term3 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * np.log(freq_ratio) * exp_arg_deriv / sc_time_freq - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 return np.sum(deriv_mod_int, axis = 2) @@ -2284,23 +2609,24 @@ def deriv2_model_wrt_scattering_index_dm_index(parameters: dict, model: float, c spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq - freq_diff = np.log(model.freqs[freq_idx]) * model.freqs[freq_idx] ** dm_index - \ + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = np.log(model.freqs[freq]) * model.freqs[freq] ** dm_index - \ np.log(ref_freq) * ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = -(1 + burst_width ** 2 / sc_time_freq ** 2 - current_timediff[freq_idx, :] / sc_time_freq) * np.log(freq_ratio) + term0 = -(1 + burst_width ** 2 / sc_time_freq ** 2 - current_timediff[freq, :] / sc_time_freq) * np.log(freq_ratio) term0_deriv = -np.log(freq_ratio) * dm_const * dm * freq_diff / sc_time_freq - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) erf_arg_deriv = -dm_const * dm * freq_diff / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) + spectrum *= freq_ratio ** (-sc_index) exp_arg_deriv = dm_const * dm * freq_diff / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = term0_deriv * model.spectrum_per_component[freq_idx, :, current_component] - term2 = term0 * deriv_first[freq_idx, :] + term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] + term2 = term0 * deriv_first[freq, :] term3 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * np.log(freq_ratio) * exp_arg_deriv / sc_time_freq - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 return np.sum(deriv_mod_int, axis = 2) @@ -2343,22 +2669,23 @@ def deriv2_model_wrt_dm_index_dm_index(parameters: dict, model: float, component spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - for freq_idx in range(num_freq): - freq_ratio = model.freqs[freq_idx] / ref_freq + for freq in range(num_freq): + freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index - erf_arg = (current_timediff[freq_idx, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq_idx, :] / sc_time_freq - erf_arg ** 2 + erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) + exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - product = dm_const * dm * (np.log(model.freqs[freq_idx]) * model.freqs[freq_idx] ** dm_index -\ + spectrum *= freq_ratio ** (-sc_index) + product = dm_const * dm * (np.log(model.freqs[freq]) * model.freqs[freq] ** dm_index -\ np.log(ref_freq) * ref_freq ** dm_index) - product_deriv = dm_const * dm * (np.log(model.freqs[freq_idx]) ** 2 * model.freqs[freq_idx] ** dm_index -\ + product_deriv = dm_const * dm * (np.log(model.freqs[freq]) ** 2 * model.freqs[freq] ** dm_index -\ np.log(ref_freq) ** 2 * ref_freq ** dm_index) exp_arg_deriv = product / sc_time_freq - 2 * erf_arg * erf_arg_deriv - term1 = product_deriv * model.spectrum_per_component[freq_idx, :, current_component] / sc_time_freq - term2 = product * deriv_first[freq_idx, :] / sc_time_freq + term1 = product_deriv * model.spectrum_per_component[freq, :, current_component] / sc_time_freq + term2 = product * deriv_first[freq, :] / sc_time_freq term3 = -np.sqrt(2 / np.pi) * product_deriv * spectrum * np.exp(exp_arg) / burst_width / sc_time_freq term4 = -np.sqrt(2 / np.pi) * product * spectrum * np.exp(exp_arg) / burst_width / sc_time_freq * exp_arg_deriv - deriv_mod_int[freq_idx, :, current_component] = term1 + term2 + term3 + term4 + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4 return np.sum(deriv_mod_int, axis = 2) From 196ba84fcfa3a6633ad4e449ff5499a1abfb0e8f Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 11 Oct 2023 09:45:28 -0700 Subject: [PATCH 5/8] initial structure for downsampling, but must be developed further. --- fitburst/pipelines/fitburst_example_chimefrb.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fitburst/pipelines/fitburst_example_chimefrb.py b/fitburst/pipelines/fitburst_example_chimefrb.py index d29544b..f11357f 100644 --- a/fitburst/pipelines/fitburst_example_chimefrb.py +++ b/fitburst/pipelines/fitburst_example_chimefrb.py @@ -65,7 +65,7 @@ ) parser.add_argument( - "--downsample_freq", action="store", dest="factor_freq_downsample", default=64, type=int, + "--downsample_freq", action="store", dest="factor_freq_downsample", default=1, type=int, help="Downsample the raw spectrum along the frequency axis by a specified integer." ) @@ -341,7 +341,7 @@ ) # if desired, downsample data prior to extraction. - data.downsample(factor_freq_downsample, factor_time_upsample) + #data.downsample(factor_freq_downsample, factor_time_upsample) log.info(f"downsampled raw data by factors of (ds_freq, ds_time) = ({factor_freq_downsample}, {factor_time_downsample})") # if the number of RFI-flagged channels is "too large", skip this event altogether. From 61790a370328390022bdd3bef5c293ddf2442d28 Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 18 Oct 2023 09:11:43 -0700 Subject: [PATCH 6/8] fixed bugs in compute_hessian, added factor of 0.5 in covariance def --- fitburst/analysis/fitter.py | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/fitburst/analysis/fitter.py b/fitburst/analysis/fitter.py index b0e9ee8..0369d74 100644 --- a/fitburst/analysis/fitter.py +++ b/fitburst/analysis/fitter.py @@ -101,28 +101,30 @@ def compute_hessian(self, data: float, parameter_list: list) -> float: # define the scale of the Hessian matrix and its labels. par_labels_output = [] + par_labels_idxs = [] par_labels = [] num_par = 0 - print("here are the fit parameters:", self.fit_parameters) for current_par in self.fit_parameters: if np.any([current_par == x for x in self.global_parameters]): par_labels_output += [current_par] + par_labels_idxs += [0] par_labels += [current_par] num_par += 1 else: par_labels_output += [f"{current_par}{idx+1}" for idx in range(self.model.num_components)] + par_labels_idxs += [idx for idx in range(self.model.num_components)] par_labels += ([current_par] * self.model.num_components) num_par += (self.model.num_components) + # now loop over all fit parameters and compute derivatives. hessian = np.zeros((num_par, num_par), dtype=float) - # now loop over all fit parameters and compute derivatives. for current_par_idx_1, current_par_1 in zip(range(num_par), par_labels): current_par_deriv_1 = getattr(deriv, f"deriv_model_wrt_{current_par_1}") current_deriv_1 = current_par_deriv_1( - parameter_dict, self.model, component=(current_par_idx_1 % self.model.num_components) + parameter_dict, self.model, component=par_labels_idxs[current_par_idx_1] ) # for efficient calculation, only compute one half of the matrix and fill the other half appropriately. @@ -131,7 +133,7 @@ def compute_hessian(self, data: float, parameter_list: list) -> float: # also compute a given derivative for all burst components. current_par_deriv_2 = getattr(deriv, f"deriv_model_wrt_{current_par_2}") current_deriv_2 = current_par_deriv_2( - parameter_dict, self.model, component=(current_par_idx_2 % self.model.num_components) + parameter_dict, self.model, component=par_labels_idxs[current_par_idx_2] ) # correct name ordering of mixed partial derivative, if necessary. @@ -141,17 +143,28 @@ def compute_hessian(self, data: float, parameter_list: list) -> float: except AttributeError: current_mixed_deriv = getattr(deriv, f"deriv2_model_wrt_{current_par_2}_{current_par_1}") - # only computed mixed derivative for parameters that describe the same component. + # only computed mixed derivative for *non-global* parameters that describe the same component, + # or any mixture of global/non-global or global-global pairs. current_deriv_mixed = 0 - if (current_par_idx_1 % self.model.num_components) == (current_par_idx_2 % self.model.num_components): + if (current_par_1 not in self.global_parameters) and (current_par_2 not in self.global_parameters) and \ + (par_labels_idxs[current_par_idx_1] != par_labels_idxs[current_par_idx_2]): + pass + + else: + # suss out correct pulse component and take mixed-partial derivative. + component = par_labels_idxs[current_par_idx_2] + + if current_par_2 in self.global_parameters: + component = par_labels_idxs[current_par_idx_1] + current_deriv_mixed = current_mixed_deriv( - parameter_dict, self.model, component=(current_par_idx_2 % self.model.num_components) + parameter_dict, self.model, component=component ) # finally, compute the hessian here. - current_hes = 2 * current_deriv_1 * current_deriv_2 - residual * current_deriv_mixed - hessian[current_par_idx_1, current_par_idx_2] = np.sum(current_hes * self.weights[:, None]) + current_hes = 2 * (current_deriv_1 * current_deriv_2 - residual * current_deriv_mixed) + hessian[current_par_idx_1, current_par_idx_2] = np.sum(current_hes * self.weights[:, None] ** 2) hessian[current_par_idx_2, current_par_idx_1] = hessian[current_par_idx_1, current_par_idx_2] return hessian, par_labels_output @@ -460,7 +473,7 @@ def _compute_fit_statistics(self, spectrum_observed: float, fit_result: object) hessian_approx = fit_result.jac.T.dot(fit_result.jac) covariance_approx = np.linalg.inv(hessian_approx) * chisq_final_reduced hessian, par_labels = self.compute_hessian(self.data, self.fit_parameters) - covariance = np.linalg.inv(hessian) * chisq_final_reduced + covariance = np.linalg.inv(0.5 * hessian) * chisq_final_reduced uncertainties = [float(x) for x in np.sqrt(np.diag(covariance)).tolist()] self.covariance_approx = covariance_approx From 92adde00d29b65c3a4c49c26465770877a6f6c91 Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 18 Oct 2023 09:12:23 -0700 Subject: [PATCH 7/8] fixed small bug in freq input to compute_profile --- fitburst/analysis/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fitburst/analysis/model.py b/fitburst/analysis/model.py index 38a226d..039a164 100644 --- a/fitburst/analysis/model.py +++ b/fitburst/analysis/model.py @@ -219,7 +219,7 @@ def compute_model(self, data: float = None) -> float: current_sc_time, current_sc_idx, current_width, - current_freq_arr, + current_freq_arr[:, None], current_ref_freq, is_folded = self.is_folded, ) From 248e7dcd42444ca5875ed19cbde6fe0d76c0a7de Mon Sep 17 00:00:00 2001 From: Emmanuel Fonseca Date: Wed, 18 Oct 2023 09:13:44 -0700 Subject: [PATCH 8/8] first set of major updates to derivative functions --- fitburst/routines/derivative.py | 231 +++++++++++++++++++------------- 1 file changed, 139 insertions(+), 92 deletions(-) diff --git a/fitburst/routines/derivative.py b/fitburst/routines/derivative.py index e14a293..2022ccc 100644 --- a/fitburst/routines/derivative.py +++ b/fitburst/routines/derivative.py @@ -272,22 +272,43 @@ def deriv2_argument_erf(name1: str, name2: str, freq: float, time_diff: float, p sc_time_freq = sc_time * (freq / ref_freq) ** sc_index deriv_mixed = 0. - if name1 == "arrival_time" and name2 == "burst_width": + if name1 == "arrival_time" and name2 == "arrival_time": + pass # mixed partial evaluates to 0. + + elif name1 == "arrival_time" and name2 == "burst_width": deriv_mixed = 1 / width ** 2 / np.sqrt(2) + elif name1 == "arrival_time" and name2 == "scattering_timescale": + pass # mixed partial evaluates to 0. + + elif name1 == "arrival_time" and name2 == "dm": + pass # mixed partial evaluates to 0. + elif name1 == "burst_width" and name2 == "burst_width": deriv_mixed = np.sqrt(2) * time_diff / width ** 3 - elif (name1 == "dm" and name2 == "burst_width") or (name1 == "dm_index" and name2 == "burst_width"): - deriv_mixed = -deriv_time_dm(name1, freq, ref_freq, dm, dm_index) / width ** 2 / np.sqrt(2) + elif (name1 == "burst_width" and name2 == "dm") or (name1 == "burst_width" and name2 == "dm_index"): + deriv_mixed = -deriv_time_dm(name2, freq, ref_freq, dm, dm_index) / width ** 2 / np.sqrt(2) - elif name1 == "scattering_index" and name2 == "burst_width": + elif name1 == "burst_width" and name2 == "scattering_index": log_freq = np.log(freq / ref_freq) deriv_mixed = log_freq / np.sqrt(2) / sc_time_freq - elif name1 == "scattering_timescale" and name2 == "burst_width": + elif name1 == "burst_width" and name2 == "scattering_timescale": deriv_mixed = 1 / sc_time / sc_time_freq / np.sqrt(2) + elif name1 == "dm" and name2 == "dm": + pass # mixed partial evaluates to 0. + + elif name1 == "dm" and name2 == "scattering_timescale": + pass # mixed partial evaluates to 0. + + elif name1 == "scattering_timescale" and name2 == "scattering_timescale": + deriv_mixed = -np.sqrt(2) * width / sc_time_freq / sc_time ** 2 + + else: + sys.exit(f"ERROR: unrecognized deriv2_argument_erf pair: {name1}, {name2}") + return deriv_mixed def deriv_model_wrt_amplitude(parameters: dict, model: object, component: int = 0) -> float: @@ -624,8 +645,7 @@ def deriv_model_wrt_dm_index(parameters: dict, model: float, component: int = 0, term1 = -deriv_timediff_wrt_dm[freq] * model.spectrum_per_component[freq, :, current_component] term1 /= sc_times_freq[freq] term2 = spectrum * np.exp(arg_exp) * deriv_arg_erf * 2 / np.sqrt(np.pi) - - deriv_mod_int[freq, :, current_component] = (term1 + term2) + deriv_mod_int[freq, :, current_component] = (term1 + term2) # now determine if all components should be summed or not. deriv_mod = deriv_mod_int[:, :, component] @@ -1593,7 +1613,7 @@ def deriv2_model_wrt_burst_width_scattering_timescale(parameters: dict, model: f "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, component ) deriv2_arg_erf = deriv2_argument_erf( - "scattering_timescale", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + "burst_width", "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, component ) # now define terms that contriubte to mixed-partial derivative. @@ -1657,7 +1677,7 @@ def deriv2_model_wrt_burst_width_scattering_index(parameters: dict, model: float "scattering_index", model.freqs[freq], current_timediff[freq, :], parameters, component ) deriv2_arg_erf = deriv2_argument_erf( - "scattering_index", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + "burst_width", "scattering_index", model.freqs[freq], current_timediff[freq, :], parameters, component ) # now define terms that contriubte to mixed-partial derivative. @@ -1715,7 +1735,7 @@ def deriv2_model_wrt_burst_width_dm(parameters: dict, model: float, component: i # if scattering is not resolvable, then assume Gaussian temporal profile. if sc_time_freq < np.fabs(0.15 * burst_width): deriv_time_dm_freq = deriv_time_dm("dm", model.freqs[freq], ref_freq, dm, dm_index) - deriv_mod[freq, :] += -2 * current_timediff[freq, :] * deriv_time_dm_freq * \ + deriv_mod[freq, :] += 2 * current_timediff[freq, :] * deriv_time_dm_freq * \ current_model[freq, :] / burst_width ** 3 deriv_mod[freq, :] += (current_timediff[freq, :] ** 2 * deriv_first[freq, :] / burst_width ** 3) @@ -1731,7 +1751,7 @@ def deriv2_model_wrt_burst_width_dm(parameters: dict, model: float, component: i "dm", model.freqs[freq], current_timediff[freq, :], parameters, component ) deriv2_arg_erf = deriv2_argument_erf( - "dm", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + "burst_width", "dm", model.freqs[freq], current_timediff[freq, :], parameters, component ) # now define terms that contriubte to mixed-partial derivative. @@ -1802,7 +1822,7 @@ def deriv2_model_wrt_burst_width_dm_index(parameters: dict, model: float, compon "dm_index", model.freqs[freq], current_timediff[freq, :], parameters, component ) deriv2_arg_erf = deriv2_argument_erf( - "dm_index", "burst_width", model.freqs[freq], current_timediff[freq, :], parameters, component + "burst_width", "dm_index", model.freqs[freq], current_timediff[freq, :], parameters, component ) # now define terms that contriubte to mixed-partial derivative. @@ -1859,14 +1879,21 @@ def deriv2_model_wrt_arrival_time_arrival_time(parameters: dict, model: float, c else: spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) spectrum *= freq_ratio ** (-sc_index) - erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) - erf_arg_deriv = -1 / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = (1 / sc_time_freq - 2 * erf_arg * erf_arg_deriv) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "arrival_time", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "arrival_time", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "arrival_time", "arrival_time", model.freqs[freq], current_timediff[freq, :], parameters, component + ) term1 = deriv_first[freq, :] / sc_time_freq - term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) - deriv_mod[freq, :] = term1 + term2 + term2 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term3 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod[freq, :] = term1 + term2 + term3 return deriv_mod @@ -1897,7 +1924,7 @@ def deriv2_model_wrt_arrival_time_dm(parameters: dict, model: float, component: current_timediff = model.timediff_per_component[:, :, component] deriv_first = deriv_model_wrt_dm(parameters, model, component, add_all = False) deriv_mod = np.zeros(current_model.shape, dtype=float) - dm_const = 4149.377593360996 + dm = parameters["dm"][0] # global parameter. dm_index = parameters["dm_index"][0] # global parameter. ref_freq = parameters["ref_freq"][component] sc_index = parameters["scattering_index"][0] @@ -1906,32 +1933,37 @@ def deriv2_model_wrt_arrival_time_dm(parameters: dict, model: float, component: spectral_running = parameters["spectral_running"][component] # now loop over each frequency and compute mixed-derivative array per channel. + deriv_timediff_wrt_dm = deriv_time_dm("dm", model.freqs, ref_freq, dm, dm_index) + for freq in range(current_model.shape[0]): freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index if sc_time_freq < np.fabs(0.15 * burst_width): deriv_mod[freq, :] = current_timediff[freq, :] * deriv_first[freq, :] / burst_width ** 2 - deriv_mod[freq, :] -= (dm_const * (model.freqs[freq] ** dm_index - ref_freq ** dm_index) * - current_model[freq, :] / burst_width ** 2) + deriv_mod[freq, :] += deriv_timediff_wrt_dm[freq] * current_model[freq, :] / burst_width ** 2 else: # adjust time-difference values to make them friendly for error function. current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) spectrum *= freq_ratio ** (-sc_index) - erf_arg = (current_timediff[freq, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = -dm_const * (model.freqs[freq] ** dm_index - ref_freq ** dm_index) / burst_width / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = dm_const * (model.freqs[freq] ** dm_index - ref_freq ** dm_index) / sc_time_freq - \ - 2 * erf_arg * erf_arg_deriv + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "arrival_time", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "dm", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "arrival_time", "dm", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. term1 = deriv_first[freq, :] / sc_time_freq - term2 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) - deriv_mod[freq, :] = term1 + term2 + term2 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term3 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod[freq, :] = term1 + term2 + term3 return deriv_mod @@ -1998,7 +2030,6 @@ def deriv2_model_wrt_arrival_time_scattering_timescale(parameters: dict, model: deriv_first = deriv_model_wrt_scattering_timescale(parameters, model, component, add_all = False) deriv_mod = np.zeros(current_model.shape) dm = parameters["dm"][0] # global parameter. - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] # global parameter. ref_freq = parameters["ref_freq"][component] dm_index = parameters["dm_index"][0] # global parameter. @@ -2016,19 +2047,22 @@ def deriv2_model_wrt_arrival_time_scattering_timescale(parameters: dict, model: sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) spectrum *= freq_ratio ** (-sc_index) - erf_arg = (current_timediff[freq, :] - burst_width ** 2 / sc_time_freq) / burst_width / np.sqrt(2) - erf_arg_deriv = burst_width * freq_ratio ** sc_index / sc_time_freq ** 2 / np.sqrt(2) - exp_arg = burst_width ** 2 / 2 / sc_time_freq ** 2 - \ - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = -burst_width ** 2 * freq_ratio ** sc_index / sc_time_freq ** 3 + \ - current_timediff[freq, :] * freq_ratio ** sc_index / sc_time_freq ** 2 - \ - 2 * erf_arg * erf_arg_deriv + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "arrival_time", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv_arg_exp = deriv_argument_exp( + "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, component + ) + deriv2_arg_erf = deriv2_argument_erf( + "arrival_time", "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, component + ) # now define terms that contriubte to mixed-partial derivative. term1 = -model.spectrum_per_component[freq, :, component] / sc_time / sc_time_freq term2 = deriv_first[freq, :] / sc_time_freq - term3 = 0. # np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * freq_ratio ** sc_index / burst_width / sc_time_freq ** 2 - term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + term3 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term4 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) deriv_mod[freq, :] = term1 + term2 + term3 + term4 return deriv_mod @@ -2061,7 +2095,6 @@ def deriv2_model_wrt_arrival_time_scattering_index(parameters: dict, model: floa deriv_first = deriv_model_wrt_scattering_index(parameters, model, component, add_all = False) deriv_mod = np.zeros(current_model.shape) dm = parameters["dm"][0] # global parameter. - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] # global parameter. ref_freq = parameters["ref_freq"][component] dm_index = parameters["dm_index"][0] # global parameter. @@ -2090,8 +2123,8 @@ def deriv2_model_wrt_arrival_time_scattering_index(parameters: dict, model: floa # now define terms that contriubte to mixed-partial derivative. term1 = -model.spectrum_per_component[freq, :, component] * np.log(freq_ratio) / sc_time_freq term2 = deriv_first[freq, :] / sc_time_freq - term3 = np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * np.log(freq_ratio) / burst_width * freq_ratio ** (-sc_index) - term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + term3 = np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * np.log(freq_ratio) / burst_width + term4 = -np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width deriv_mod[freq, :] = term1 + term2 + term3 + term4 return deriv_mod @@ -2156,7 +2189,7 @@ def deriv2_model_wrt_dm_dm(parameters: dict, model: float, component: int = 0) - """ deriv_mod_int = np.zeros(model.spectrum_per_component.shape, dtype=float) - dm_const = 4149.377593360996 + dm = parameters["dm"][0] # global parameter. dm_index = parameters["dm_index"][0] # global parameter. sc_index = parameters["scattering_index"][0] sc_time = parameters["scattering_timescale"][0] @@ -2171,32 +2204,39 @@ def deriv2_model_wrt_dm_dm(parameters: dict, model: float, component: int = 0) - spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] - # now loop over each frequency and compute mixed-derivative array per channel. - for freq in range(current_model.shape[0]): - freq_ratio = model.freqs[freq] / ref_freq - freq_diff = model.freqs[freq] ** dm_index - ref_freq ** dm_index - sc_time_freq = sc_time * freq_ratio ** sc_index + # now loop over each frequency and compute mixed-derivative array per channel. + deriv_timediff_wrt_dm = deriv_time_dm("dm", model.freqs, ref_freq, dm, dm_index) - if sc_time_freq < np.fabs(0.15 * burst_width): - deriv_mod_int[freq, :, current_component] = (dm_const * freq_diff) ** 2 * current_timediff[freq, :] ** 2 * \ - current_model[freq, :] / burst_width ** 4 - deriv_mod_int[freq, :, current_component] -= (dm_const * freq_diff) ** 2 * current_model[freq, :] / \ - burst_width ** 2 + for freq in range(current_model.shape[0]): + freq_ratio = model.freqs[freq] / ref_freq + freq_diff = model.freqs[freq] ** dm_index - ref_freq ** dm_index + sc_time_freq = sc_time * freq_ratio ** sc_index - else: - # adjust time-difference values to make them friendly for error function. - current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width - erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) - erf_arg_deriv = -dm_const * freq_diff / burst_width / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 - spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) - spectrum *= freq_ratio ** (-sc_index) - exp_arg_deriv = dm_const * freq_diff / sc_time_freq - 2 * erf_arg * erf_arg_deriv + if sc_time_freq < np.fabs(0.15 * burst_width): + deriv_mod_int[freq, :, current_component] = (deriv_timediff_wrt_dm[freq] * current_timediff[freq, :]) ** 2 * \ + current_model[freq, :] / burst_width ** 4 + deriv_mod_int[freq, :, current_component] -= (deriv_timediff_wrt_dm[freq] / burst_width) ** 2 * current_model[freq, :] + + else: + # adjust time-difference values to make them friendly for error function. + current_timediff[current_timediff < -5 * burst_width] = -5 * burst_width + spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) + spectrum *= freq_ratio ** (-sc_index) + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "dm", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + deriv_arg_exp = deriv_argument_exp( + "dm", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + deriv2_arg_erf = deriv2_argument_erf( + "dm", "dm", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) - term1 = dm_const * freq_diff * deriv_first[freq, :] / sc_time_freq - term2 = -dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width * \ - freq_ratio ** (-sc_index) * exp_arg_deriv - deriv_mod_int[freq, :, current_component] = term1 + term2 + term1 = -deriv_timediff_wrt_dm[freq] * deriv_first[freq, :] / sc_time_freq + term2 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term3 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) + deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 return np.sum(deriv_mod_int, axis = 2) @@ -2222,7 +2262,6 @@ def deriv2_model_wrt_scattering_timescale_scattering_timescale(parameters: dict, # get dimensions and define empty model-derivative matrix. dm = parameters["dm"][0] - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] num_freq, num_time, num_component = model.spectrum_per_component.shape sc_index = parameters["scattering_index"][0] @@ -2241,23 +2280,27 @@ def deriv2_model_wrt_scattering_timescale_scattering_timescale(parameters: dict, for freq in range(num_freq): freq_ratio = model.freqs[freq] / ref_freq sc_time_freq = sc_time * freq_ratio ** sc_index - term0 = (-burst_width ** 2 / sc_time_freq ** 2 / sc_time + current_timediff[freq, :] / \ - sc_time / sc_time_freq) - term0_deriv = (3 * burst_width ** 2 / (sc_time * sc_time_freq) ** 2 - - 2 * current_timediff[freq, :] / sc_time ** 2 / sc_time_freq) - erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) - erf_arg_deriv = burst_width / sc_time_freq / sc_time / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) spectrum *= freq_ratio ** (-sc_index) - exp_arg_deriv = -(burst_width / sc_time_freq )** 2 / sc_time + current_timediff[freq, :] / sc_time_freq / sc_time - \ - 2 * erf_arg * erf_arg_deriv - + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + deriv_arg_exp = deriv_argument_exp( + "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + deriv2_arg_erf = deriv2_argument_erf( + "scattering_timescale", "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + # now compute derivative terms here. + term0 = (current_timediff[freq, :] / sc_time_freq - (burst_width / sc_time_freq) ** 2) / sc_time + term0_deriv = -(current_timediff[freq, :] / sc_time_freq - (burst_width / sc_time_freq) ** 2) / sc_time ** 2 + \ + (-current_timediff[freq, :] / sc_time_freq + 2 * (burst_width / sc_time_freq) ** 2) / sc_time ** 2 term1 = term0_deriv * model.spectrum_per_component[freq, :, current_component] term2 = term0 * deriv_first[freq, :] - term3 = 0. # -np.sqrt(8 / np.pi) * spectrum * burst_width * np.exp(exp_arg) / sc_time_freq / sc_time ** 2 - term4 = np.sqrt(2 / np.pi) * spectrum * burst_width * np.exp(exp_arg) * exp_arg_deriv / sc_time_freq + term3 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term4 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4 @@ -2285,7 +2328,6 @@ def deriv2_model_wrt_scattering_timescale_scattering_index(parameters: dict, mod # get dimensions and define empty model-derivative matrix. dm = parameters["dm"][0] - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] num_freq, num_time, num_component = model.spectrum_per_component.shape sc_index = parameters["scattering_index"][0] @@ -2348,7 +2390,6 @@ def deriv2_model_wrt_scattering_timescale_dm(parameters: dict, model: float, com # get dimensions and define empty model-derivative matrix. dm = parameters["dm"][0] - dm_const = 4149.377593360996 dm_index = parameters["dm_index"][0] num_freq, num_time, num_component = model.spectrum_per_component.shape sc_index = parameters["scattering_index"][0] @@ -2364,23 +2405,29 @@ def deriv2_model_wrt_scattering_timescale_dm(parameters: dict, model: float, com ref_freq = parameters["ref_freq"][current_component] spectral_index = parameters["spectral_index"][current_component] spectral_running = parameters["spectral_running"][current_component] + deriv_timediff_wrt_dm = deriv_time_dm("dm", model.freqs, ref_freq, dm, dm_index) + # TODO: replace lines below. for freq in range(num_freq): freq_ratio = model.freqs[freq] / ref_freq - freq_diff = model.freqs[freq] ** dm_index - ref_freq ** dm_index sc_time_freq = sc_time * freq_ratio ** sc_index spectrum = 10 ** amplitude * freq_ratio ** (spectral_index + spectral_running * np.log(freq_ratio)) spectrum *= freq_ratio ** (-sc_index) - erf_arg = (current_timediff[freq, :] / burst_width - burst_width / sc_time_freq) / np.sqrt(2) - erf_arg_deriv = burst_width / sc_time / sc_time_freq / np.sqrt(2) - exp_arg = (burst_width / sc_time_freq) ** 2 / 2 - current_timediff[freq, :] / sc_time_freq - erf_arg ** 2 - exp_arg_deriv = -1 / sc_time * (burst_width / sc_time_freq) ** 2 + current_timediff[freq, :] / sc_time / sc_time_freq \ - -2 * erf_arg * erf_arg_deriv + arg_exp = argument_exp(burst_width, current_timediff[freq, :], sc_time_freq) + deriv_arg_erf = deriv_argument_erf( + "dm", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + deriv_arg_exp = deriv_argument_exp( + "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) + deriv2_arg_erf = deriv2_argument_erf( + "dm", "scattering_timescale", model.freqs[freq], current_timediff[freq, :], parameters, current_component + ) - term1 = dm_const * freq_diff * deriv_first[freq, :] / sc_time_freq - term2 = -dm_const * freq_diff * current_model[freq, :] / sc_time / sc_time_freq - term3 = dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) / burst_width / freq_ratio ** (-sc_index) / sc_time - term4 = -dm_const * freq_diff * np.sqrt(2 / np.pi) * spectrum * np.exp(exp_arg) * exp_arg_deriv / burst_width * freq_ratio ** (-sc_index) + term1 = -deriv_timediff_wrt_dm[freq] * deriv_first[freq, :] / sc_time_freq + term2 = deriv_timediff_wrt_dm[freq] * current_model[freq, :] / sc_time / sc_time_freq + term3 = spectrum * np.exp(arg_exp) * deriv2_arg_erf * 2 / np.sqrt(np.pi) + term4 = spectrum * np.exp(arg_exp) * deriv_arg_erf * deriv_arg_exp * 2 / np.sqrt(np.pi) deriv_mod_int[freq, :, current_component] = term1 + term2 + term3 + term4