Commit 8052f691 authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Refactored decomposition code

parent 4c654b37

etl/assets/data/source_properties.yml

0 → 100644
+28 −0
Changes for etl/assets/data/source_properties.yml: 28 added lines, 0 removed lines.
Original line number Diff line number Diff line 
G332.96: {

  CH3OH: {

    'data': 'G332.96_CH3OH_3mm_feather.fits',

    'rms': 0.0034,

    'vlsr': -67.6,

    'frequency_limits_ghz': [ 96.755, 96.775 ] }

}

G28: {

  CH3OH: {

    'data': 'G28_CH3OH_3mm_feather.fits',

    'rms': 0.003,

    'vlsr': 79.6,

    'frequency_limits_ghz': [ 96.6995, 96.740 ] }

}

G14: {

  CH3OH: {

    'data': 'G14_CH3OH_3mm_feather.fits',

    'rms': 0.0032,

    'vlsr': 40.5,

    'frequency_limits_ghz': [ 96.706, 96.750 ] }

}

G12: {

  CH3OH: {

    'data': 'G12_CH3OH_3mm_ACA_MA.fits',

    'rms': 0.002,

    'vlsr': 34.4,

    'frequency_limits_ghz': [ 96.692, 96.786 ] }

}
 
 No newline at end of file

etl/prs/prs_visfit.py

+85 −59
Changes for etl/prs/prs_visfit.py: 85 added lines, 59 removed lines.
Original line number Diff line number Diff line 
import pickle

import plotly.graph_objects as go

import plotly.express as px

import plotly.offline as poff

import matplotlib.pyplot as plt

import numpy as np





source_names = ['G332', 'G28', 'G14', 'G12']

source = source_names[0]

def convert_ratio_to_density(ratio):

    return 10 ** ((np.arctanh((ratio - 0.523) / 1.009) + 6.086) ** 0.971)



with open(f'final_results_{source}_fullres.pickle', 'rb') as infile:

    fit_results = pickle.load(infile)



rms = 0.003

detection_86 = fit_results[0,:,:] > 3 * rms

detection_A = fit_results[3,:,:] > 3 * rms

discard_failed_86 = (fit_results[0,:,:] < 0.1) & (fit_results[2,:,:] > 0.3)

discard_failed_87 = fit_results[4,:,:] < 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



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)

ncomps = np.nansum(~np.isnan(fit_results[(0, 6, 12),:,:]), axis=0)

peak_component = np.nanmax(fit_results[(0, 6, 12),:,:], axis=0)

comp_idx = np.argmax(fit_results[(0, 6, 12),:,:] == peak_component, axis=0)

velocity_field = np.where((comp_idx == 0) & (fit_results[1,:,:] > -200), fit_results[1,:,:], velocity_field_1)

velocity_field = np.where((comp_idx == 1) & (fit_results[7,:,:] > -200), fit_results[7,:,:], velocity_field)

velocity_field = np.where((comp_idx == 2) & (fit_results[13,:,:] > -200), fit_results[13,:,:], velocity_field)

velocity_field = np.where(velocity_field > -200, velocity_field, np.nan)



grid = np.meshgrid(np.arange(512),np.arange(512))

fig = go.Figure(data=[go.Scatter3d(x=grid[0].flatten(), y=grid[1].flatten(),z=velocity_field_1.flatten(), mode='markers', opacity=0.1, marker=dict(size=1,symbol='cross'))])

fig.add_scatter3d(x=grid[0].flatten(), y=grid[1].flatten(),z=velocity_field_2.flatten(), mode='markers', opacity=0.1, marker=dict(size=1,symbol='cross'))

fig.add_scatter3d(x=grid[0].flatten(), y=grid[1].flatten(),z=velocity_field_3.flatten(), mode='markers', opacity=0.1, marker=dict(size=1,symbol='cross'))



fig = go.Figure(data=[go.Surface(z=velocity_field)])



ratio_1 = np.where(detection_86 & discard_failed & unsaturated_87_86,

                   fit_results[4,:,:] / fit_results[0,:,:],

                   np.nan)

ratio_2 = np.where((fit_results[6,:,:] > 3 * 0.003) & (fit_results[6,:,:] < 1) & (fit_results[10,:,:] < 1) & (fit_results[10,:,:] / fit_results[6,:,:] < 1.1), fit_results[10,:,:] / fit_results[6,:,:], np.nan)

ratio_3 = np.where((fit_results[12,:,:] > 3 * 0.003) & (fit_results[12,:,:] < 1) & (fit_results[16,:,:] < 1) & (fit_results[16,:,:] / fit_results[12,:,:] < 1.1), fit_results[16,:,:] / fit_results[12,:,:], np.nan)

ratio_map = np.where((comp_idx == 0) & detection_86 & discard_failed & unsaturated_87_86, fit_results[4,:,:] / fit_results[0,:,:], np.nan)

ratio_map = np.where((comp_idx == 1) & (fit_results[6,:,:] > 3 * 0.003) & discard_failed & (fit_results[10,:,:] / fit_results[6,:,:] < 1.2), fit_results[10,:,:] / fit_results[6,:,:], ratio_map)

ratio_map = np.where((comp_idx == 2) & (fit_results[12,:,:] > 3 * 0.003) & discard_failed & (fit_results[16,:,:] / fit_results[12,:,:] < 1.2), fit_results[16,:,:] / fit_results[12,:,:], ratio_map)

ratio_map = np.where(fit_results[0,:,:] < fit_results[6,:,:], ratio_2, ratio_1)

density_map = 10 ** ((np.arctanh(ratio_map - 0.55) + 2.3) ** 0.9 + 3.75)

density_map = np.where(density_map < 5e4, 5e4, density_map)

def rotate_z(x, y, z, theta):

    w = x + 1j * y

    return np.real(np.exp(1j * theta) * w), np.imag(np.exp(1j * theta) * w), z



size_arr = np.where(~np.isnan(density_map), fit_results[0,:,:], 0)

fig = go.Figure(data=[go.Surface(z=density_map, surfacecolor=size_arr)])

# grid = np.meshgrid(np.arange(512),np.arange(512))

# fig = go.Figure(data=[go.Scatter3d(x=grid[0].flatten(), y=grid[1].flatten(),z=density_map.flatten(), mode='markers', opacity=0.3, marker=dict(size=size_arr, symbol='cross'))])



x_eye = -1.25

y_eye = 2

z_eye = 2

def make_rotating_figure(fig,

                         marker_size=2,

                         marker_color='blue',

                         marker_opacity=0.5,

                         x_eye=-1.25,

                         y_eye=2,

                         z_eye=2,

                         outfile=None):

    _outfile = outfile if outfile else 'rotating_3d_scatter.html'

    fig.update_traces(marker=dict(size=marker_size, color=marker_color, opacity=marker_opacity, symbol='circle'))

    fig.update_layout(

        title=None,

        width=600,
@@ -64,7 +32,6 @@ fig.update_layout( 
        paper_bgcolor='rgba(0,0,0,0)',

        plot_bgcolor='rgba(0,0,0,0)',

        font_color="white",

         scene=dict(zaxis=dict(type='log', dtick=1)),

        margin=dict(l=0, r=0, b=0, t=0),

        modebar=dict(bgcolor='rgba(0, 0, 0, 0)'),

        scene_camera_eye=dict(x=x_eye, y=y_eye, z=z_eye),
@@ -77,7 +44,7 @@ fig.update_layout( 
                          pad=dict(t=0, r=0),

                          buttons=[dict(label='Play',

                                        method='animate',

                                visible=False,

                                        visible=True,

                                        args=[None, dict(frame=dict(duration=5, redraw=True),

                                                         transition=dict(duration=0),

                                                         fromcurrent=True,
@@ -88,17 +55,76 @@ fig.update_layout( 
                          )

                     ]

    )





def rotate_z(x, y, z, theta):

    w = x+1j*y

    return np.real(np.exp(1j*theta)*w), np.imag(np.exp(1j*theta)*w), z



    frames = []

    for t in np.arange(0, 6.26, 0.1):

        xe, ye, ze = rotate_z(x_eye, y_eye, z_eye, -t)

        frames.append(go.Frame(layout=dict(scene_camera_eye=dict(x=xe, y=ye, z=ze))))

    fig.frames = frames

plotly.offline.plot(fig, filename=r'test_animation.html')

    poff.plot(fig, filename=_outfile)

    # fig.show()





source_names = ['G332.96', 'G28', 'G14', 'G12']

source_names = ['G28',]



rms = 0.003

marker_opacity = 0.05

marker_size = 1.5

marker_color = 'DarkSlateGrey'

species = 'CH3OH'



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)

    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]))



    fig = px.scatter_3d()

    for velocity in [velocity_field_1, velocity_field_2, velocity_field_3]:

        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,

                       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,

                       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,

                       fit_results[16, :, :] / fit_results[12, :, :], np.nan)

    fig = px.scatter_3d(log_z=True)

    for ratio in [ratio_1, ratio_2, ratio_3]:

        mask = ~np.isnan(ratio.flatten())

        fig.add_scatter3d(x=grid[0].flatten()[mask], y=grid[1].flatten()[mask], z=convert_ratio_to_density(ratio).flatten()[mask],

                          mode='markers')



    make_rotating_figure(fig=fig, z_eye=1, outfile=f'rotating_3d_scatter_density_{source}_{species}.html',

                         marker_size=marker_size, marker_color=marker_color, marker_opacity=marker_opacity)

fig.show()

etl/requirements.txt

+1 −0
Changes for etl/requirements.txt: 1 added line, 0 removed lines.
Original line number Diff line number Diff line
@@ -7,3 +7,4 @@ astroquery 
pandas

pyspeckit

spectral-cube

plotly
 
 No newline at end of file