Revision of decomposition script to unify CH3OH and N2H+ (d45528d5) · Commits · Andrea Giannetti / CoCoA_methanol_archive

etl/mdl/mdl_psk_decomposition.py

+44 −39

Changes for etl/mdl/mdl_psk_decomposition.py: 44 added lines, 39 removed lines.
Original line number	Diff line number	Diff line
		@@ -29,7 +29,7 @@ def get_bic(residuals: np.array,
		return bic


		def convert_dvelocity_to_hz(dvelocity):
		def convert_dvelocity_to_hz(dvelocity, rest_frequencies):
		return dvelocity / speed_of_light.to(u.km / u.s) * rest_frequencies[0].to(u.Hz)


		@@ -39,10 +39,10 @@ def create_limits(par_number, channelwidth_hz, rms):
		if idx % 3 == 0:
		limits.append((rms, 10))
		elif idx % 3 == 2:
		limits.append((channelwidth_hz, 20))
		limits.append((channelwidth_hz/2, 7.5e5))
		else:
		limits.append((0, 0))
		limited = [(True, False) if (idx % 3 == 0) or (idx % 3 == 2) else (False, False) for idx in range(par_number)]
		limited = [(True, True) if (idx % 3 == 0) or (idx % 3 == 2) else (False, False) for idx in range(par_number)]
		return limits, limited


		@@ -52,31 +52,30 @@ def fill_tied_pars(tied_pars):
		return _tied_pars_list


		def fill_guesses(line_peak, doppler_adjusted_frequencies, shift_hz, linewidth_guess_hz, rms):
		guesses = [line_peak, doppler_adjusted_frequencies[0].to(u.Hz).value + shift_hz, linewidth_guess_hz,
		line_peak, doppler_adjusted_frequencies[1].to(u.Hz).value + shift_hz, linewidth_guess_hz,
		np.max((0.3 * line_peak, rms)), doppler_adjusted_frequencies[2].to(u.Hz).value + shift_hz,
		linewidth_guess_hz,
		np.max((0.1 * line_peak, rms)), doppler_adjusted_frequencies[3].to(u.Hz).value + shift_hz,
		def fill_guesses(line_peak, doppler_adjusted_frequencies, shift_hz, linewidth_guess_hz, rms, peak_scaling_factors=(1, 1, 0.5, 0.3)):
		guesses = []
		for idx, freq in enumerate(doppler_adjusted_frequencies):
		guesses += [np.max((rms, peak_scaling_factors[idx] * line_peak)),
		freq.to(u.Hz).value + shift_hz,
		linewidth_guess_hz]
		return guesses


		def update_results(fitted_cube, final_pars, indices_to_keep, mask, rms, components, bic_ratio):
		def update_results(fitted_cube, final_pars, indices_to_keep, mask, rms, components, bic_ratio, rest_frequencies, independent_pars_per_component, pars_per_component):
		_final_pars = final_pars.copy()
		fit_results = fitted_cube.parcube.copy()
		detections = mask.copy()
		for idx_to_change in range(components + 1):
		velocity_map = (fitted_cube.parcube[1 + (12 * idx_to_change), :, :] * u.Hz).to(u.km / u.s,
		velocity_map = (fitted_cube.parcube[1 + (pars_per_component * idx_to_change), :, :] * u.Hz).to(u.km / u.s,
		equivalencies=freq_to_vel).value
		linewidth_map = 2.355 * ((fitted_cube.parcube[2 + (12 * idx_to_change), :, :] * u.Hz) / rest_frequencies[0].to(
		linewidth_map = 2.355 * ((fitted_cube.parcube[2 + (pars_per_component * idx_to_change), :, :] * u.Hz) / rest_frequencies[0].to(
		u.Hz)).value * speed_of_light.to(u.km / u.s).value
		fit_results[1 + (12 * idx_to_change), :, :] = velocity_map
		fit_results[2 + (12 * idx_to_change), :, :] = linewidth_map
		detections = (detections) & (fit_results[0 + (12 * idx_to_change), :, :] >= 2 * rms)
		_final_pars[0:6 + 6 * components, :, :] = np.where((bic_ratio > 1) & mask & detections,
		fit_results[1 + (pars_per_component * idx_to_change), :, :] = velocity_map
		fit_results[2 + (pars_per_component * idx_to_change), :, :] = linewidth_map
		detections = (detections) & (fit_results[0 + (pars_per_component * idx_to_change), :, :] >= 2 * rms)
		_final_pars[0:independent_pars_per_component + independent_pars_per_component * components, :, :] = np.where((bic_ratio > 1) & mask & detections,
		fit_results[indices_to_keep, :, :],
		final_pars[0:6 + 6 * components, :, :])
		final_pars[0:independent_pars_per_component + independent_pars_per_component * components, :, :])
		return _final_pars


		@@ -130,14 +129,14 @@ datadir = os.path.join('assets', 'data')
		with open(os.path.join(datadir, 'source_properties.yml'), 'r') as infile:
		source_properties = yaml.load(infile, Loader=yaml.FullLoader)

		species = 'CH3OH'
		species = 'N2H+'
		source_names = ['G332.96', 'G28', 'G14']
		source_names = ['G332.96']
		linewidth_guess = 1.5 * u.km / u.s
		# source_names = ['G28']
		linewidth_guess = 0.8 * u.km / u.s
		max_components = 3
		rest_frequencies = {
		'CH3OH': np.array([96.739358, 96.74142, 96.744545, 96.755501]) * u.GHz,
		'N2H+': np.array([93.1958,]) * u.GHz
		'N2H+': np.array([93.1762595,]) * u.GHz
		}

		for source in source_names:
		@@ -155,33 +154,34 @@ for source in source_names:
		smoothed_data[i, :, :] = scipy.ndimage.gaussian_filter(cube.cube[i, :, :], sigma=3)
		cube.cube = smoothed_data

		freq_to_vel = u.doppler_radio(rest_frequencies[0])
		frequency_offset = vlsr.to(u.Hz, equivalencies=freq_to_vel) - rest_frequencies[0]
		doppler_adjusted_frequencies = rest_frequencies + frequency_offset
		freq_to_vel = u.doppler_radio(rest_frequencies[species][0])
		frequency_offset = vlsr.to(u.Hz, equivalencies=freq_to_vel) - rest_frequencies[species][0]
		doppler_adjusted_frequencies = rest_frequencies[species] + frequency_offset
		if frequency_limits is not None:
		cube = cube.slice(frequency_limits[0], frequency_limits[1], unit=u.GHz)
		linewidth_guess_hz = convert_dvelocity_to_hz(dvelocity=linewidth_guess).value / 2.355
		linewidth_guess_hz = convert_dvelocity_to_hz(dvelocity=linewidth_guess,
		rest_frequencies=rest_frequencies[species]).value / 2.355
		channelwidth_hz = cube.header['CDELT3']
		if len(rest_frequencies) > 1:
		delta_freq = rest_frequencies[1:] - rest_frequencies[0]
		if len(rest_frequencies[species]) > 1:
		delta_freq = rest_frequencies[species][1:] - rest_frequencies[species][0]
		independent_indices = np.array((0, 1, 2, 3, 6, 9))
		indices_to_keep = np.array([], dtype=int)
		independent_pars_per_component = len(independent_indices)
		channel_separation = delta_freq.to(u.Hz).value / np.abs(channelwidth_hz)
		else:
		delta_freq = rest_frequencies[1:] - rest_frequencies[0]
		independent_indices = np.array((0, 1, 2, 3, 6, 9))
		independent_indices = np.array((0, 1, 2,))
		indices_to_keep = np.array([], dtype=int)
		independent_pars_per_component = len(independent_indices)
		pars_per_component = 3 * len(rest_frequencies)
		pars_per_component = 3 * len(rest_frequencies[species])
		final_pars = np.zeros((pars_per_component * max_components, cube.cube.shape[1], cube.cube.shape[2])) * np.nan



		threshold_spec_sum = 4 * np.sqrt(cube.shape[0]) * rms
		if 'starting_point' in source_properties[source][species]:
		starting_point = source_properties[source][species]['starting_point']
		else:
		spec_sum = np.nansum(cube.cube, axis=0)
		# spec_sum = -np.where(spec_sum < 0, 0, spec_sum)
		# spec_sum[230:240, 295:305] *= -1
		threshold_spec_sum = 4 * np.sqrt(cube.shape[0]) * rms
		starting_point = np.unravel_index(np.argmax(spec_sum), spec_sum.shape)
		mask = np.where(spec_sum > threshold_spec_sum, True, False)
		bic = -1e10 * np.ones_like(spec_sum)
		@@ -196,7 +196,7 @@ for source in source_names:
		if np.sum(mask) > 0:
		idx_shift = components * pars_per_component

		if len(rest_frequencies) > 1:
		if len(rest_frequencies[species]) > 1:
		tied_pars = {
		'PAR_VELO': 1 + idx_shift,
		'DELTA_FREQ_1': delta_freq[0].to(u.Hz).value,
		@@ -206,7 +206,7 @@ for source in source_names:
		}
		tied += fill_tied_pars(tied_pars=tied_pars)
		else:
		tied = None
		tied = [''] * pars_per_component
		# best_guesses = get_guesses(peaks=peaks,
		# ea_channel_separation=channel_separation[0],
		# ref_freq=cube.header['CRVAL3'],
		@@ -219,7 +219,7 @@ for source in source_names:
		shift_hz=linewidth_guess_hz / 2 * idx_shift,
		linewidth_guess_hz=linewidth_guess_hz,
		rms=rms)
		limits, limited = create_limits(par_number=len(tied),
		limits, limited = create_limits(par_number=len(guesses),
		channelwidth_hz=np.abs(channelwidth_hz),
		rms=rms)
		indices_to_keep = np.concatenate((indices_to_keep, (independent_indices + idx_shift)))
		@@ -232,7 +232,7 @@ for source in source_names:
		limits=limits,
		limited=limited,
		maskmap=mask,
		use_nearest_as_guess=True,
		# use_nearest_as_guess=True,
		use_neighbor_as_guess=True,
		multicore=10)
		# cube.plot_spectrum(x=starting_point[1], y=starting_point[0])
		@@ -254,10 +254,15 @@ for source in source_names:
		mask=mask,
		rms=rms,
		components=components,
		bic_ratio=bic / new_bic)
		bic_ratio=bic / new_bic,
		independent_pars_per_component=independent_pars_per_component,
		pars_per_component=pars_per_component,
		rest_frequencies=rest_frequencies[species])
		print(bic / new_bic)
		bic = new_bic
		threshold = (1.1 * cube.shape[0] + cube.shape[0] * np.log(2 * np.pi)
		threshold = (3 * (np.sqrt(2 * np.abs(linewidth_guess_hz / channelwidth_hz) * (components + 1) * len(rest_frequencies[species])))
		+ cube.shape[0]
		+ cube.shape[0] * np.log(2 * np.pi)
		+ cube.shape[0] * np.log(rms ** 2) + np.log(cube.shape[0]) * number_of_parameters)
		mask = np.where(bic > threshold, True, False)

etl/prs/prs_visfit.py

+24 −28

Changes for etl/prs/prs_visfit.py: 24 added lines, 28 removed lines.
Original line number	Diff line number	Diff line
		@@ -4,6 +4,8 @@ import plotly.express as px
		import plotly.offline as poff
		import matplotlib.pyplot as plt
		import numpy as np
		import matplotlib
		matplotlib.use('TkAgg')


		def convert_ratio_to_density(ratio):
		@@ -71,53 +73,47 @@ rms = 0.003
		marker_opacity = 0.05
		marker_size = 1.5
		marker_color = 'DarkSlateGrey'
		species = 'CH3OH'
		species = 'N2H+'

		main_line_intensity_indices = (3, 9, 15) if species == 'CH3OH' else (0, 3, 6)
		velocity_indices = (1, 7, 13) if species == 'CH3OH' else (1, 4, 7)

		for source in source_names:
		with open(f'final_results_{source}_{species}.pickle', 'rb') as infile:
		fit_results = pickle.load(infile)

		components_map = np.nansum(fit_results[(3, 9, 15), :, :] > 3 * rms, axis=0)
		components_map = np.nansum(fit_results[main_line_intensity_indices, :, :] > 3 * rms, axis=0)
		plt.clf()
		plt.imshow(components_map, origin='lower')
		plt.colorbar()
		detection_86 = fit_results[0, :, :] > 3 * rms
		detection_86_2comps = fit_results[6, :, :] > 3 * rms
		detection_86_3comps = fit_results[12, :, :] > 3 * rms
		detection_A = fit_results[3, :, :] > 3 * rms
		discard_failed_86 = (fit_results[0, :, :] < 0.1) & (fit_results[2, :, :] > 0.3)
		discard_failed_86_2comps = (fit_results[6, :, :] < 0.1) & (fit_results[2, :, :] > 0.3)
		discard_failed_86_3comps = (fit_results[12, :, :] < 0.1) & (fit_results[2, :, :] > 0.3)
		discard_failed_87 = fit_results[4, :, :] < 0.1
		discard_failed_87_2comps = fit_results[10, :, :] < 0.1
		discard_failed_87_3comps = fit_results[16, :, :] < 0.1
		discard_failed_88 = fit_results[5, :, :] < 0.1
		discard_failed = discard_failed_86 & discard_failed_87 & discard_failed_88
		unsaturated_87_86 = fit_results[4, :, :] / fit_results[0, :, :] < 1.2
		unsaturated_87_86_2comps = fit_results[10, :, :] / fit_results[6, :, :] < 1.2
		unsaturated_87_86_3comps = fit_results[16, :, :] / fit_results[12, :, :] < 1.2
		velocity_field_1 = np.where(detection_A & discard_failed_86 & (np.abs(fit_results[1, :, :]) < 200),
		fit_results[1, :, :], np.nan)
		velocity_field_2 = np.where(detection_A & discard_failed_86 & (np.abs(fit_results[7, :, :]) < 200),
		fit_results[7, :, :], np.nan)
		velocity_field_3 = np.where(detection_A & discard_failed_86 & (np.abs(fit_results[13, :, :]) < 200),
		fit_results[13, :, :], np.nan)
		grid = np.meshgrid(np.arange(detection_A.shape[0]), np.arange(detection_A.shape[0]))
		detection_main_line = fit_results[main_line_intensity_indices[0], :, :] > 3 * rms
		velocity_field = []
		for component in range(3):
		velocity_field.append(np.where(detection_main_line, fit_results[velocity_indices[component], :, :], np.nan))
		grid = np.meshgrid(np.arange(detection_main_line.shape[0]), np.arange(detection_main_line.shape[0]))

		fig = px.scatter_3d()
		for velocity in [velocity_field_1, velocity_field_2, velocity_field_3]:
		for velocity in velocity_field:
		mask = ~np.isnan(velocity.flatten())
		fig.add_scatter3d(x=grid[0].flatten()[mask], y=grid[1].flatten()[mask], z=velocity.flatten()[mask],
		mode='markers')
		make_rotating_figure(fig=fig, z_eye=0.5, outfile=f'rotating_3d_scatter_velocity_{source}_{species}.html',
		marker_size=marker_size, marker_color=marker_color, marker_opacity=marker_opacity)

		ratio_1 = np.where(detection_86 & discard_failed & unsaturated_87_86,
		if species == 'CH3OH':
		detection_86 = fit_results[0, :, :] > 3 * rms
		detection_87 = (fit_results[4, :, :] > rms) & (fit_results[4, :, :] * fit_results[2, :, :] * 2.355 > 3 * (rms * np.sqrt(2 * fit_results[2, :, :] * 2.355 / 0.3)))
		detection_86_2comps = fit_results[6, :, :] > 3 * rms
		detection_86_3comps = fit_results[12, :, :] > 3 * rms
		unsaturated_87_86 = fit_results[4, :, :] / fit_results[0, :, :] < 1.2
		unsaturated_87_86_2comps = fit_results[10, :, :] / fit_results[6, :, :] < 1.2
		unsaturated_87_86_3comps = fit_results[16, :, :] / fit_results[12, :, :] < 1.2
		ratio_1 = np.where(detection_86 & detection_87 & unsaturated_87_86,
		fit_results[4, :, :] / fit_results[0, :, :],
		np.nan)
		ratio_2 = np.where(detection_86_2comps & discard_failed_86_2comps & discard_failed_87_2comps & unsaturated_87_86_2comps,
		ratio_2 = np.where(detection_86_2comps & unsaturated_87_86_2comps,
		fit_results[10, :, :] / fit_results[6, :, :], np.nan)
		ratio_3 = np.where(detection_86_3comps & discard_failed_86_3comps & discard_failed_87_3comps & unsaturated_87_86_3comps,
		ratio_3 = np.where(detection_86_3comps & unsaturated_87_86_3comps,
		fit_results[16, :, :] / fit_results[12, :, :], np.nan)
		fig = px.scatter_3d(log_z=True)
		for ratio in [ratio_1, ratio_2, ratio_3]: