<jupyter_start><jupyter_text>Method 1: Using PystanI'll do this for a single n, and treat the mode locations as fixed for now.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import lightkurve as lk
from astropy.units import cds
from astropy import units as u
import corner
import matplotlib
import pystan
import pickle
import os
matplotlib.rc('xtick', labelsize=20)
matplotlib.rc('ytick', labelsize=20)
matplotlib.rc('axes',labelsize=20)
matplotlib.rcParams['text.usetex'] = False

#The stellar parameters are for 16 Cyg A
nus = 0.411 #uHz
i = np.deg2rad(56.) #rad
d02 = 6.8 #uHz
dnu = 102. #uHz
locs = np.genfromtxt('locs.txt')
mid = int(np.floor(len(locs)/2))
l0 = locs[mid, 0]
l1 = locs[mid, 1]
l2 = locs[mid-1, 2]
l3 = locs[mid-1, 3]
ff = np.genfromtxt('freqs.txt')
pp = np.genfromtxt('model.txt')
sel = [(ff >= 2185) & (ff <= 2200)]
f = ff[tuple(sel)]
p = pp[tuple(sel)] 
# modelocs = np.concatenate((l0, l1, l2, l3))
modelocs = np.array([l2, l0])
# modeids = [0]*len(l0) + [1]*len(l1) + [2]*len(l2) + [3]*len(l3)
modeids = [2, 0]
pg = lk.Periodogram(f*u.microhertz, p*(cds.ppm**2/u.microhertz))
ax = pg.plot(alpha=.5)
pg.smooth(filter_width=2.).plot(ax=ax, linewidth=2)
plt.scatter(modelocs, [15]*len(modelocs),c=modeids, s=20, edgecolor='k')
plt.show()<jupyter_output>/home/oliver/.local/lib/python3.6/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))<jupyter_text>Lets fit a group of lorentzians to these data!<jupyter_code>overwrite = True
pbjam = '''
functions{      
    vector lorentzian(real loc, int l, int m, vector f, real eps, real H, real w, real nus){
        return (eps * H) ./ (1 + (4/w^2) * square(f - loc + m*nus));
    }
}
data{
    int N;            // Number of data points
    int M;            // Number of modes
    vector[N] f;      // Frequency
    vector[N] p;      // Power
    real asy_locs[M]; // Mode locations (this will have to change for multiple n modes)
    int asy_ids[M];   // The ID's of the modes
}
parameters{
    real logAmp[M];         // Mode amplitude in log space
    real logGamma[M];       // Mode linewidth in log space
    real locs[M];  // True mode locations
    real vsini;    // Line of sight rotational frequency
    real<lower=0.> nus;      // Rotational frequency splitting
}
transformed parameters{
    real H[M];       // Mode height
    real w[M];       // Mode linewidth
    
    for (m in 1:M){
        w[m] = 10^logGamma[m];             // Transform mode linewidth from log space
        H[m] = 10^logAmp[m] / pi() / w[m]; // Transform mode amplitude to mode height
    }
}
model{
    vector<lower=1.>[N] modes; // Our Model
    real i;          // Angle of inclination (rad)
    matrix[4,4] eps; // Matrix of legendre polynomials
    int l;           // The radial degree
    
    // First we'll calculate all the legendre polynomials for this i
    i = asin(vsini / nus);
    eps = rep_matrix(1., 4, 4);
    eps[0+1,0+1] = 1.;
    eps[1+1,0+1] = cos(i)^2;
    eps[1+1,1+1] = 0.5 * sin(i)^2;
    eps[2+1,0+1] = 0.25 * (3. * cos(i)^2 - 1.)^2;
    eps[2+1,1+1] = (3./8.)*sin(2*i)^2;
    eps[2+1,2+1] = (3./8.) * sin(i)^4;
    eps[3+1,0+1] = (1./64.)*(5.*cos(3.*i) + 3.*cos(i))^2;
    eps[3+1,1+1] = (3./64.)*(5.*cos(2.*i) + 3.)^2 * sin(i)^2;
    eps[3+1,2+1] = (15./8.)*cos(i)^2 * sin(i)^4;
    eps[3+1,3+1] = (5./16.)*sin(i)^6;    

    modes = rep_vector(1., N);
    for (mode in 1:M){        // Iterate over all modes passed in
        l = asy_ids[mode];    // Identify the Mode ID
        for (m in -l:l){      // Iterate over all m in a given l
            modes += lorentzian(locs[mode], l, m, f, eps[l+1,abs(m)+1], H[mode], w[mode], nus);
        }
    }

    // Model drawn from a gamma distribution scaled to the model (Anderson+1990)
    p ~ gamma(1., 1../modes);
    
    //priors on the parameters
    logAmp ~ normal(1.5, 1);
    logGamma ~ normal(0, 0.01);
    locs ~ normal(asy_locs, 1);
    nus ~ normal(0.411, 0.1);
    vsini ~ uniform(0.,nus);
}
'''
model_path = 'pbjam.pkl'
if overwrite:
    print('Updating Stan model')
    sm = pystan.StanModel(model_code = pbjam, model_name='pbjam')
    pkl_file =  open(model_path, 'wb')
    pickle.dump(sm, pkl_file)
    pkl_file.close()
if os.path.isfile(model_path):
    print('Reading in Stan model')
    sm = pickle.load(open(model_path, 'rb'))
else:
    print('Saving Stan Model')
    sm = pystan.StanModel(model_code = pbjam, model_name='pbjam')
    pkl_file =  open(model_path, 'wb')
    pickle.dump(sm, pkl_file)
    pkl_file.close()
data = {'N':len(f),
        'M': len(modelocs),
        'f':f,
        'p':p,
        'asy_locs':modelocs,
        'asy_ids':modeids}

init = {'logAmp' :   np.ones(len(modelocs))*1.5,
        'logGamma' : np.zeros(len(modelocs)),
        'nus':nus,
        'locs' : modelocs}
nchains = 4
fit = sm.sampling(data=data, iter=50, chains=nchains, init = [init for n in range(nchains)])
print(fit)
sini = fit['vsini']/fit['nus']
import seaborn as sns
sns.distplot(sini)
test = fit['test']
test.shape
e10 = test[:,1, 0]
e11 = test[:,1, 1]
e20 = test[:,2, 0]
e21 = test[:,2, 1]
e22 = test[:,2, 2]
e30 = test[:,3, 0]
e31 = test[:,3, 1]
e32 = test[:,3, 2]
e33 = test[:,3, 3]
labels = ['10','11','20','21','22','30','31','32','33']
chains = np.array([e10,e11,e20,e21,e22,e30,e31,e32,e33])
import corner
corner.corner(chains.T,labels=labels, quantiles=[0.16, 0.5, 0.84],show_titles=True)
plt.show()<jupyter_output>WARNING:root:Too few points to create valid contours
WARNING:root:Too few points to create valid contours
{'fraction_non_alphanumeric': 0.17934589406327764, 'fraction_numerical': 0.049768929968005686, 'mean_word_length': 3.0953420669577874, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '563691', 'license': 'MIT', 'n_tokens_mistral': 2427, 'n_tokens_neox': 2255, 'n_words': 556, 'path': 'code/tests/completed_tests/test_i_causes_nans.ipynb', 'repo_name': 'ojhall94/malatium', 'size': '5626'}