Skip to content

Commit

Permalink
allow for all integers (including 0) in fourth column as instrument ID;
Browse files Browse the repository at this point in the history 
fixes #63; fix some instrument names in plots
  • Loading branch information
@j-faria
j-faria committed Jun 11, 2020
1 parent c27b190 commit 1316a1a
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 22 deletions.
35 changes: 21 additions & 14 deletions pykima/display.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1072,7 +1072,7 @@ def kima_pars_to_keplerian_pars(p):
**ekwargs)
ax.fill_between(t[m], -jitters[k - 1], jitters[k - 1],
alpha=0.2)
print(get_instrument_name(self.data_file[k-1]), end=': ')
print(self.instruments[k-1], end=': ')
print(wrms(self.y[m] - v[m] - KOvel[m] - GPvel[m], 1/e[m]**2))
residuals[m] = self.y[m] - v[m] - KOvel[m] - GPvel[m]

Expand Down Expand Up @@ -1446,8 +1446,7 @@ def make_plot5(self, show=True, ranges=None):
units = ['m/s'] * self.n_jitters + ['m/s', 'days', 'days', None]
xlabels = []
for i in range(self.n_jitters):
l = r'%s$_{{\rm %s}}}$' % (labels[i],
get_instrument_name(self.data_file[i]))
l = r'%s$_{{\rm %s}}}$' % (labels[i], self.instruments[i])

xlabels.append(l)
for label, unit in zip(labels[self.n_jitters:],
Expand Down Expand Up @@ -1904,13 +1903,14 @@ def plot_random_planets(self, ncurves=50, over=0.1, pmin=None, pmax=None,
else:
if self.multi:
for j in range(self.inst_offsets.shape[1] + 1):
inst = self.instruments[j]
m = self.obs == j + 1
ax.errorbar(t[m], y[m], yerr[m], fmt='o', color=colors[j],
label=os.path.basename(self.data_file[j]))
label=inst)
if self.KO:
ax1.errorbar(t[m], y[m] - v_KO_at_t.mean(axis=0)[m],
yerr[m], fmt='o', color=colors[j],
label=os.path.basename(self.data_file[j]))
label=inst)
ax.legend(loc='upper left')
else:
ax.errorbar(t, y, yerr, fmt='o')
Expand Down Expand Up @@ -2278,8 +2278,10 @@ def hist_vsys(self, show_offsets=True, specific=None):
axs = [axs,]

for i in range(n_inst_offsets):
wrt = get_instrument_name(self.data_file[-1])
this = get_instrument_name(self.data_file[i])
# wrt = get_instrument_name(self.data_file[-1])
# this = get_instrument_name(self.data_file[i])
wrt = self.instruments[-1]
this = self.instruments[i]
label = 'offset\n%s rel. to %s' % (this, wrt)
a = self.inst_offsets[:, i]
estimate = percentile68_ranges_latex(a) + units
Expand All @@ -2306,17 +2308,21 @@ def hist_vsys(self, show_offsets=True, specific=None):
i = specific.index(self.data_file[-1])
# toggle: if i is 0 it becomes 1, if it's 1 it becomes 0
i ^= 1
wrt = get_instrument_name(self.data_file[-1])
this = get_instrument_name(specific[i])
# wrt = get_instrument_name(self.data_file[-1])
# this = get_instrument_name(specific[i])
wrt = self.instruments[-1]
this = self.instruments[i]
label = 'offset\n%s rel. to %s' % (this, wrt)
offset = self.inst_offsets[:, self.data_file.index(specific[i])]
estimate = percentile68_ranges_latex(offset) + units
fig, ax = plt.subplots(1, 1, constrained_layout=True)
ax.hist(offset)
ax.set(xlabel=label, title=estimate, ylabel='posterior samples')
else:
wrt = get_instrument_name(specific[1])
this = get_instrument_name(specific[0])
# wrt = get_instrument_name(specific[1])
# this = get_instrument_name(specific[0])
wrt = self.instruments[specific[1]]
this = self.instruments[specific[0]]
label = 'offset\n%s rel. to %s' % (this, wrt)
of1 = self.inst_offsets[:, self.data_file.index(specific[0])]
of2 = self.inst_offsets[:, self.data_file.index(specific[1])]
Expand All @@ -2339,10 +2345,11 @@ def hist_extra_sigma(self):
if self.multi: # there are n_instruments jitters
# lambda substrs
fig, axs = plt.subplots(1, self.n_instruments, sharey=True,
figsize=(self.n_instruments * 3,
5), squeeze=True)
figsize=(self.n_instruments * 3, 5),
squeeze=True)
for i, jit in enumerate(self.extra_sigma.T):
inst = get_instrument_name(self.data_file[i])
# inst = get_instrument_name(self.data_file[i])
inst = self.instruments[i]
estimate = percentile68_ranges_latex(jit) + units
axs[i].hist(jit)
axs[i].set(xlabel='%s' % inst, title=estimate,
Expand Down
8 changes: 8 additions & 0 deletions pykima/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,14 @@ def read_datafile(datafile, skip):
else:
data = np.loadtxt(datafile, usecols=(0, 1, 2), skiprows=skip)
obs = np.loadtxt(datafile, usecols=(3, ), skiprows=skip, dtype=int)
uobs = np.unique(obs)

id0 = 0
for i, o in enumerate(obs):
if o != uobs[id0]:
id0 += 1
obs[i] = id0 + 1

return data, obs


Expand Down
33 changes: 25 additions & 8 deletions src/Data.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,14 +248,31 @@ void Data::load_multi(const std::string filename, const std::string units, int s
double factor = 1.;
if(units == "kms") factor = 1E3;

for (unsigned n = 0; n < data.size(); n++)
{
if (n<skip) continue;
t.push_back(data[n][0]);
y.push_back(data[n][1] * factor);
sig.push_back(data[n][2] * factor);
obsi.push_back(data[n][3]);
}
// the 4th column of the file identifies the instrument; can have "0"s
// this is to make sure the obsi vector always starts at 1, to avoid
// segmentation faults later
vector<int> inst_id;
int id = 0;

inst_id.push_back(data[skip][3]);

for (size_t n = skip+1; n < data.size(); n++)
{
if (data[n][3] != inst_id.back())
inst_id.push_back(data[n][3]);
}

for (unsigned n = skip; n < data.size(); n++)
{
// if (n < skip) continue;
t.push_back(data[n][0]);
y.push_back(data[n][1] * factor);
sig.push_back(data[n][2] * factor);

if (data[n][3] != inst_id[id])
id++;
obsi.push_back(id+1);
}

// How many points did we read?
printf("# Loaded %zu data points from file %s\n", t.size(), filename.c_str());
Expand Down

0 comments on commit 1316a1a

Please sign in to comment.