<commit_before>import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import emcee
from matplotlib.widgets import Slider, Button, RadioButtons
from scipy.stats import spearmanr,pearsonr,kendalltau
from scipy.optimize import minimize
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d,UnivariateSpline
import pdb
import pandas as pd
import time
from scipy import signal

def air_to_vacuum(airwl,nouvconv=True):
    """
    Returns vacuum wavelength of the provided air wavelength array or scalar.
    Good to ~ .0005 angstroms.

    If nouvconv is True, does nothing for air wavelength < 2000 angstroms.
    
    Input must be in angstroms.
    
    Adapted from idlutils airtovac.pro, based on the IAU standard 
    for conversion in Morton (1991 Ap.J. Suppl. 77, 119)
    """
    airwl = np.array(airwl,copy=False,dtype=float,ndmin=1)
    isscal = airwl.shape == tuple()
    if isscal:
        airwl = airwl.ravel()
    
    #wavenumber squared
    sig2 = (1e4/airwl)**2
    
    convfact = 1. + 6.4328e-5 + 2.94981e-2/(146. - sig2) +  2.5540e-4/( 41. - sig2)
    newwl = airwl.copy() 
    if nouvconv:
        convmask = newwl>=2000
        newwl[convmask] *= convfact[convmask]
    else:
        newwl[:] *= convfact
    return newwl[0] if isscal else newwl

def gaussian_lines(line_x,line_a,xgrid,width=2.0):
    '''
    Creates ideal Xenon spectrum
    '''
    #print 'Creating ideal calibration spectrum'
    temp = np.zeros(xgrid.size)
    for i in range(line_a.size):
        gauss = line_a[i]*np.exp(-(xgrid-line_x[i])**2/(2*width**2))
        temp += gauss
    return temp

def wavecalibrate(px,fx,slit_x,stretch_est=None,shift_est=None,qu_es=None):
    #flip and normalize flux
    fx = fx - np.min(fx)
    fx = fx[::-1]
    fx = fx/signal.medfilt(fx,201)

    #prep calibration lines into 1d spectra
    wm,fm = np.loadtxt('osmos_Xenon.dat',usecols=(0,2),unpack=True)
    wm = air_to_vacuum(wm)
    xgrid = np.arange(0.0,6800.0,0.01)
    lines_gauss = gaussian_lines(wm,fm,xgrid)
    interp = interp1d(xgrid,lines_gauss,bounds_error=False,fill_value=0)
    #interp = UnivariateSpline(xgrid,lines_gauss)

    wave_est = qu_es*(px-slit_x)**2+(px)*stretch_est+shift_est
    wm_in = wm[np.where((wm<wave_est.max())&(wm>wave_est.min()))]
    px_max = np.zeros(wm_in.size)
    for i in range(wm_in.size):
        px_in = px[np.where((wave_est<wm_in[i]+5.0)&(wave_est>wm_in[i]-5))]
        px_max[i] = px_in[fx[np.where((wave_est<wm_in[i]+5.0)&(wave_est>wm_in[i]-5))].argmax()]

    def polyfour(x,a,b,c,d,e,f):
        return a + b*x + c*x**2.0 + d*x**3.0 + e*x**4.0 + f*x**5.0

    params,pcov = curve_fit(polyfour,(px_max-slit_x),wm_in,p0=[2000.0,0.70,1e-5,1e-8,1e-12,1e-12])
    '''
    def log_prior(theta):
        #g_i needs to be between 0 and 1
        if (all(theta > 0) and all(theta < 1)):
            return 0
        else:
            return -np.inf  # recall log(0) = -inf

    def log_likelihood(p, x, y, e):
        dy = y - (p[0]+p[1]*x + p[2]*x*x + p[3]*x*x*x + p[4]*x*x*x*x + p[5]*x*x*x*x*x)
        #g = np.clip(p[6:], 0, 1)  # g<0 or g>1 leads to NaNs in logarithm
        logL1 = - 0.5 * np.log(2 * np.pi * e ** 2) - 0.5 * (dy / e) ** 2<commit_msg>changed to fitting within 5000Ang<commit_after>import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import emcee
from matplotlib.widgets import Slider, Button, RadioButtons
from scipy.stats import spearmanr,pearsonr,kendalltau
from scipy.optimize import minimize
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d,UnivariateSpline
import pdb
import pandas as pd
import time
from scipy import signal

def air_to_vacuum(airwl,nouvconv=True):
    """
    Returns vacuum wavelength of the provided air wavelength array or scalar.
    Good to ~ .0005 angstroms.

    If nouvconv is True, does nothing for air wavelength < 2000 angstroms.
    
    Input must be in angstroms.
    
    Adapted from idlutils airtovac.pro, based on the IAU standard 
    for conversion in Morton (1991 Ap.J. Suppl. 77, 119)
    """
    airwl = np.array(airwl,copy=False,dtype=float,ndmin=1)
    isscal = airwl.shape == tuple()
    if isscal:
        airwl = airwl.ravel()
    
    #wavenumber squared
    sig2 = (1e4/airwl)**2
    
    convfact = 1. + 6.4328e-5 + 2.94981e-2/(146. - sig2) +  2.5540e-4/( 41. - sig2)
    newwl = airwl.copy() 
    if nouvconv:
        convmask = newwl>=2000
        newwl[convmask] *= convfact[convmask]
    else:
        newwl[:] *= convfact
    return newwl[0] if isscal else newwl

def gaussian_lines(line_x,line_a,xgrid,width=2.0):
    '''
    Creates ideal Xenon spectrum
    '''
    #print 'Creating ideal calibration spectrum'
    temp = np.zeros(xgrid.size)
    for i in range(line_a.size):
        gauss = line_a[i]*np.exp(-(xgrid-line_x[i])**2/(2*width**2))
        temp += gauss
    return temp

def wavecalibrate(px,fx,slit_x,stretch_est=None,shift_est=None,qu_es=None):
    #flip and normalize flux
    fx = fx - np.min(fx)
    fx = fx[::-1]
    fx = fx/signal.medfilt(fx,201)

    #prep calibration lines into 1d spectra
    wm,fm = np.loadtxt('osmos_Xenon.dat',usecols=(0,2),unpack=True)
    wm = air_to_vacuum(wm)
    xgrid = np.arange(0.0,6800.0,0.01)
    lines_gauss = gaussian_lines(wm,fm,xgrid)
    interp = interp1d(xgrid,lines_gauss,bounds_error=False,fill_value=0)
    #interp = UnivariateSpline(xgrid,lines_gauss)

    wave_est = qu_es*(px-slit_x)**2+(px)*stretch_est+shift_est
    wm_in = wm[np.where((wm<wave_est.max())&(wm>wave_est.min()))]
    wm_in = wm[np.where((wm<5000.0)&(wm>wave_est.min()))]
    px_max = np.zeros(wm_in.size)
    for i in range(wm_in.size):
        px_in = px[np.where((wave_est<wm_in[i]+5.0)&(wave_est>wm_in[i]-5))]
        px_max[i] = px_in[fx[np.where((wave_est<wm_in[i]+5.0)&(wave_est>wm_in[i]-5))].argmax()]

    def polyfour(x,a,b,c,d,e,f):
        return a + b*x + c*x**2.0 + d*x**3.0 + e*x**4.0 + f*x**5.0

    params,pcov = curve_fit(polyfour,(px_max-slit_x),wm_in,p0=[2000.0,0.70,1e-5,1e-8,1e-12,1e-12])
    '''
    def log_prior(theta):
        #g_i needs to be between 0 and 1
        if (all(theta > 0) and all(theta < 1)):
            return 0
        else:
            return -np.inf  # recall log(0) = -inf

    def log_likelihood(p, x, y, e):
        dy = y - (p[0]+p[1]*x + p[2]*x*x + p[3]*x*x*x + p[4]*x*x*x*x + p[5]*x*x*x*x*x)
        #g = np.clip(p[6:], 0, 1)  # g<0 or g>1 leads to NaNs in logarithm
        logL1 = - 0.5 * np.log(2 * np.pi * e ** 2) - 0.5 * (dy / e) ** 2
{'fraction_non_alphanumeric': 0.1485496183206107, 'fraction_numerical': 0.04381679389312977, 'mean_word_length': 3.1123666038920277, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, '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}}
{'dir': 'git-commits-cleaned', 'id': '1585191', 'n_tokens_mistral': 2883, 'n_tokens_neox': 2696, 'n_words': 553}