<jupyter_start><jupyter_text>Find bimodal samplings:<jupyter_code>from hq.samples_analysis import is_n_modal
bimodal_mask = np.zeros(len(meta), dtype=bool)
bimodal_modes = {}
bimodal_mode_fractions = {}
for n, row in enumerate(tqdm(meta)):
    if row['unimodal']:
        continue
        
    samples = load_samples(c, row)
    data = get_rvdata(allvisit[allvisit['APOGEE_ID'] == row['APOGEE_ID']])
    bimodal, modes, n_per_mode = is_n_modal(data, samples, n_clusters=2)
    bimodal_mask[n] = bimodal
    
    if bimodal:
        bimodal_modes[row['APOGEE_ID']] = modes
        bimodal_mode_fractions[row['APOGEE_ID']] = np.array(n_per_mode) / sum(n_per_mode)
bimodal_mask.sum()
rows = []
units = dict()
for apid in tqdm(bimodal_modes.keys()):
    for samples, frac in zip(bimodal_modes[apid],
                             bimodal_mode_fractions[apid]):
        row = {'APOGEE_ID': apid}
        
        if len(samples) == 1:
            MAP_sample = samples
        else:
            MAP_idx = (samples['ln_prior'] + samples['ln_likelihood']).argmax()
            MAP_sample = samples[MAP_idx]
            
        for k in MAP_sample.par_names:
            row['MAP_' + k] = MAP_sample[k].value
            units['MAP_' + k] = MAP_sample[k].unit
        row['t0_bmjd'] = MAP_sample.t0.tcb.mjd
        row['bimodal_mode_fraction'] = frac
        
        rows.append(row)
        
tbl = at.QTable(rows)
for k in units:
    tbl[k] = tbl[k] * units[k]
bimodal = at.join(tbl, meta, keys='APOGEE_ID', 
                  uniq_col_name='{table_name}{col_name}', 
                  table_names=['', 'meta_'])

# recompute m2_mins
percentiles = [1, 5, 16, 50, 84, 95, 99]
for row in tqdm(bimodal):
    sample = extract_MAP_sample(row)
    mf = fast_mf(sample['P'], sample['K'], sample['e'])
    m1s = np.random.normal(row['mass'], row['mass_err'], size=1000)
    m2_mins = fast_m2_min(m1s, mf.value)
    for pp in percentiles:
        row[f'm2_min_{pp}'] = np.percentile(m2_mins, pp) * u.Msun

bimodal.write('../catalogs/bimodal.fits', overwrite=True)<jupyter_output><empty_output><jupyter_text>---<jupyter_code>llr = meta['max_unmarginalized_ln_likelihood'] - meta['robust_constant_ln_likelihood']
gold = meta[~np.isin(meta['APOGEE_ID'].astype(str), reject_ids)
            & (meta['LOGG'] > -0.5) & (meta['LOGG'] < 5.5)
            & (meta['TEFF'] > 3500) & (meta['TEFF'] < 8000)
            & (meta['M_H'] > -2.5) & (meta['M_H'] < 0.5)
            & (meta['ALPHA_M'] > -2.) & (meta['ALPHA_M'] < 2.)
            & (meta['MAP_s'] < 0.5)
            & (meta['n_visits'] > 5)
            & (llr > 25)
            # & ((meta['ASPCAPFLAG'] & 2**7) == 0)
            # & (meta['phase_coverage'] > 0.5)
            # & (meta['periods_spanned'] > 1.)
            # & (meta['max_phase_gap'] < 0.5)]
           ]
gold = gold[gold['mcmc_success']]

len(gold)
gold.write('../catalogs/gold_sample.fits', overwrite=True)<jupyter_output><empty_output><jupyter_text>---<jupyter_code>plt.hist(gold['m2_min_50'] / gold['mass'], 
         bins=np.linspace(0, 1., 32));
sub_gold = gold[gold['LOGG'] > 3.5]
bin_edges = np.array([3500, 5200, 6000, 7500])
bin_ctr = 0.5 * (bin_edges[:-1] + bin_edges[1:])
bin_idx = np.digitize(sub_gold['TEFF'], bin_edges)
q_bins = np.linspace(0, 1, 16)
for i in np.arange(1, len(bin_ctr)+1):
    gold_bin = sub_gold[bin_idx == i]
    plt.figure()
    plt.hist(gold_bin['m2_min_50'] / gold_bin['mass'], 
             bins=q_bins, density=True)
# for row in gold[(gold['MAP_P'] < 10) & 
#                 (gold['LOGG'] > 4) & 
#                 (gold['MAP_e'] > 0.5)][:10]:
#     apid = row['APOGEE_ID']
#     with h5py.File(c.mcmc_results_path, 'r') as f:
#         samples = tj.JokerSamples.read(f[apid])
#         data = get_rvdata(allvisit[allvisit['APOGEE_ID'] == apid])
        
#         fig, ax = plt.subplots()
#         _ = tj.plot_rv_curves(samples[:10], data=data, ax=ax)
        
#         llr = row['max_unmarginalized_ln_likelihood'] - row['robust_constant_ln_likelihood']
#         ax.set_title(f"{row['STARFLAGS']}\n{row['ASPCAPFLAGS']}\n{row['APOGEE_ID']}: {llr:.2f}", 
#                      fontsize=12)
#         print(row['APOGEE_ID'])
# for apid in gold[(gold['MAP_P'] < 100) & (gold['LOGG'] < 1.5)]['APOGEE_ID'][:10]:
#     with h5py.File(c.mcmc_results_path, 'r') as f:
#         samples = tj.JokerSamples.read(f[apid])
#         data = get_rvdata(allvisit[allvisit['APOGEE_ID'] == apid])
        
#         fig, ax = plt.subplots()
#         _ = tj.plot_rv_curves(samples[:10], data=data, ax=ax)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1792452830188679, 'fraction_numerical': 0.026107942079859586, 'mean_word_length': 2.870118845500849, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '216274', 'license': 'MIT', 'n_tokens_mistral': 1956, 'n_tokens_neox': 1801, 'n_words': 270, 'path': 'notebooks/Make-gold-sample.ipynb', 'repo_name': 'adrn/apogee-dr16-vac', 'size': '4558'}