-
Notifications
You must be signed in to change notification settings - Fork 0
/
traingp.py
147 lines (118 loc) · 4.84 KB
/
traingp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#!/usr/bin/python
import constants
from myscience import funcs
import numpy as np
def create_train_data(velocity_axis=None,
velocity_range=None,
nchannels=None,
width_range=[5, 20],
temp_range=[40, 8000],
amp_range=None,
mean_range=None,
ncomponents=4,
nspectra=100,
rms=0.01):
'''
Parameters
----------
velocity_axis : array-like
Array of velocity value at each channel. Units of km/s.
velocity_range : array-like
Used if velocity_axis is None. Lower and upper limit to velocity range.
Paired with nchannels. Units of km/s.
nchannels : int, optional
Used if velocity_axis is None. Number of channels in spectra. Paired
with velocity_range.
width_range : array-like
Lower and upper limit to standard deviation of Gaussians expected in
spectra. Units of km/s.
temp_range : array-like
Lower and upper limit to kinetic temperatures of gas, e.g. CNM and WNM
temperatures. Temperatures converted to be width_range. Units of K.
amp_range : array-like
Lower and upper limit to amplitudes of Gaussian components in K.
mean_range : array-like
Lower and upper limit to centers of Gaussian components in km/s.
Default is (0.25 and 0.75 of (maximum velocity - minimum velocity)) +
minimum velocity.
RMS : float
RMS of spectra. Units of K / (km / s)
ncomponents : int
Number of Gaussian components in spectra.
nspectra : int
Number of spectra to generate.
Returns
-------
agd_train : dict
AGD training data set to be used as input for gausspy.gp.train function.
'''
# Get FHWM range
# --------------
if temp_range is not None:
FWHM_range = funcs.Tkin_to_FWHM(np.array(temp_range))
elif width_range is not None:
FWHM_range = funcs.std_to_FWHM(np.array(width_range))
else:
temp_range = [40, 8000]
FWHM_range = funcs.Tkin_to_FWHM(np.array(temp_range))
# Check velocity axis
# -------------------
if velocity_axis is None and velocity_range is not None:
velocity_axis = np.linspace(velocity_range[0],
velocity_range[1],
nchannels)
print 'Created velocity axis'
else:
nchannels = len(velocity_axis)
# Gaussian mean range
if mean_range is not None:
#velocity_diff = np.max(velocity_axis) - np.min(velocity_axis)
#mean_range = np.array((0.25 * velocity_diff, 0.75 * velocity_diff))
mean_range += np.min(velocity_axis)
agd_data = format_train_data(RMS=rms,
NCOMPS=ncomponents,
NCHANNELS=nchannels,
NSPECTRA=nspectra,
AMP_lims=amp_range,
FWHM_lims=FWHM_range,
MEAN_lims=mean_range,
x_value=velocity_axis)
return agd_data
def format_train_data(RMS=0.05, NCOMPS=4, NCHANNELS=512, NSPECTRA=10,
AMP_lims=None,
FWHM_lims=None,
MEAN_lims=None,
x_value=None):
# Component properties
AMP_lims = [RMS * 5, RMS * 25]
FWHM_lims = [10, 35] # channels
MEAN_lims = [0.25 * NCHANNELS, 0.75 * NCHANNELS]
if x_value is not None:
NCHANNELS = len(x_value)
# Initialize
agd_data = {}
chan = np.arange(NCHANNELS)
errors = chan * 0. + RMS # Constant noise for all spectra
# Begin populating data
for i in range(NSPECTRA):
spectrum_i = np.random.randn(NCHANNELS) * RMS
# Sample random components:
amps = np.random.rand(NCOMPS) * (AMP_lims[1] - AMP_lims[0]) + \
AMP_lims[0]
fwhms = np.random.rand(NCOMPS) * (FWHM_lims[1] - FWHM_lims[0]) + \
FWHM_lims[0]
means = np.random.rand(NCOMPS) * (MEAN_lims[1] - MEAN_lims[0]) + \
MEAN_lims[0]
# Create spectrum
for a, w, m in zip(amps, fwhms, means):
spectrum_i += funcs.gaussian(a, w, m)(chan)
# Enter results into AGD dataset
agd_data['data_list'] = agd_data.get('data_list', []) + [spectrum_i]
#agd_data['x_values'] = agd_data.get('x_values', []) + [chan]
agd_data['x_values'] = agd_data.get('x_values', []) + [x_value]
agd_data['errors'] = agd_data.get('errors', []) + [errors]
# If training data, keep answers
agd_data['amplitudes'] = agd_data.get('amplitudes', []) + [amps]
agd_data['fwhms'] = agd_data.get('fwhms', []) + [fwhms]
agd_data['means'] = agd_data.get('means', []) + [means]
return agd_data