CTN_BayesianAnalysis_Final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  # For plotting\n",
    "from matplotlib.backends.backend_pdf import PdfPages       #For saving figures to single pdf\n",
    "figlist = []\n",
    "#*******************************************************************************************************\n",
    "#Setting RC Parameters for figure size and fontsizes\n",
    "import matplotlib.pylab as pylab                           \n",
    "params = {'figure.figsize': (16, 12),\n",
    "          'xtick.labelsize':'xx-large',\n",
    "          'ytick.labelsize':'xx-large',\n",
    "          'text.usetex': False,\n",
    "          'lines.linewidth': 4,\n",
    "          'font.family': 'serif',\n",
    "          'font.serif': 'Georgia',\n",
    "          'font.size': 20,\n",
    "          'xtick.direction': 'in',\n",
    "          'ytick.direction': 'in',\n",
    "          'xtick.labelsize': 'xx-large',\n",
    "          'ytick.labelsize': 'xx-large',\n",
    "          'axes.labelsize': 'xx-large',\n",
    "          'axes.titlesize':'medium',\n",
    "          'axes.grid.axis': 'both',\n",
    "          'axes.grid.which': 'both',\n",
    "          'axes.grid': True,\n",
    "          'grid.color': 'xkcd:cement',\n",
    "          'grid.alpha': 0.3,\n",
    "          'lines.markersize': 12,\n",
    "          'lines.linewidth': 2.0,\n",
    "          'legend.borderpad': 0.2,\n",
    "          'legend.fancybox': True,\n",
    "          'legend.fontsize': 'large',\n",
    "          'legend.framealpha': 0.8,\n",
    "          'legend.handletextpad': 0.5,\n",
    "          'legend.labelspacing': 0.33,\n",
    "          'legend.loc': 'best',\n",
    "          'savefig.dpi': 140,\n",
    "          'savefig.bbox': 'tight',\n",
    "          'pdf.compression': 9}\n",
    "pylab.rcParams.update(params)\n",
    "#********************************************************************************************************\n",
    "from noiseBudgetModule import noiseBudget\n",
    "import numpy as np\n",
    "from uncertainties import ufloat as uf\n",
    "from uncertainties import unumpy as unp\n",
    "import scipy.constants as scc\n",
    "from scipy import signal\n",
    "from scipy.stats import skewnorm\n",
    "from scipy.optimize import curve_fit\n",
    "import os\n",
    "from os.path import expanduser as eu\n",
    "import time\n",
    "from collections.abc import Iterable\n",
    "from IPython.display import clear_output\n",
    "import sys\n",
    "import yaml\n",
    "import pickle\n",
    "figlist = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "\n",
    "## Prior Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normal distributed prior based on Penn et al. measurements\n",
    "def priorNphiB(phiB):\n",
    "    priorMean = 5.33e-4\n",
    "    priorStd = 2.7e-4\n",
    "    prefac = 1/np.sqrt(2*np.pi)/priorStd\n",
    "    exp = np.exp(-0.5*((priorMean -phiB)/priorStd)**2)\n",
    "    return -0.5*((priorMean -phiB)/priorStd)**2\n",
    "    # return prefac*exp\n",
    "\n",
    "def priorNphiS(phiS):\n",
    "    priorMean = 2.6e-7\n",
    "    priorStd = 2.6e-7\n",
    "    prefac = 1/np.sqrt(2*np.pi)/priorStd\n",
    "    exp = np.exp(-0.5*((priorMean -phiS)/priorStd)**2)\n",
    "    return -0.5*((priorMean -phiS)/priorStd)**2\n",
    "    # return prefac*exp\n",
    "\n",
    "# Uniformly distributed prior\n",
    "def priorUphiB(phiB):\n",
    "    if phiB>0 and phiB<200e-5:\n",
    "        P_B = 0 \n",
    "    else:\n",
    "        P_B = -np.inf\n",
    "    return P_B\n",
    "def priorUphiS(phiS):\n",
    "    if phiS>0 and phiS<94e-7:\n",
    "        P_S = 0\n",
    "    else:\n",
    "        P_S = -np.inf\n",
    "    return P_S\n",
    "\n",
    "def logprior(BulkLA, ShearLA):\n",
    "    X, Y = np.meshgrid(BulkLA, ShearLA)\n",
    "    logpriorDist = np.zeros(np.shape(X))\n",
    "    for ii in range(len(ShearLA)):\n",
    "        for jj in range(len(BulkLA)):\n",
    "            logpriorDist[ii][jj] = priorNphiB(BulkLA[jj])# + priorNphiS(ShearLA[ii])\n",
    "    return logpriorDist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Likelihood Distribution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time Series to PSD distribution conversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tsDatatoWelchArray(tsData, noverlap=None, timeSeg=5, rebinSize=5, skip=3,\n",
    "                       lowerCutoffFreq=70, upperCutoffFreq=600,\n",
    "                       average='median'):\n",
    "    timeSeries = tsData[:, 0]\n",
    "    noiseSig = tsData[:, 1]\n",
    "    SampleRate = 1/(timeSeries[1] - timeSeries[0])\n",
    "    nperseg = int(timeSeg * SampleRate)\n",
    "    #Sort data into rows with number in each row = to nperseg\n",
    "    # length of columns\n",
    "    col_len=0\n",
    "    end_edge=nperseg\n",
    "    if noverlap is None:\n",
    "        noverlap = nperseg//2\n",
    "    while(end_edge<len(noiseSig)):\n",
    "        end_edge=end_edge+nperseg-noverlap\n",
    "        col_len=col_len+1\n",
    "\n",
    "    # make array of row length nperseg and column length however much data\n",
    "    # will fit fully in row\n",
    "    sorted_data = np.zeros((col_len,nperseg))\n",
    "\n",
    "    # make array for PSD will have same col len but\n",
    "    # row len = max freq/freq_spacing + 1\n",
    "    # or (samplerate/2) /(samplerate/nperseg) + 1 which = nperseg/2 +1\n",
    "    welch_array = np.zeros((col_len, nperseg//2 + 1))\n",
    "\n",
    "    #Move through signal data array populating array\n",
    "    next_start=0\n",
    "    for ii in range(col_len):\n",
    "        sorted_data[ii, :] = noiseSig[next_start: next_start+nperseg]\n",
    "        next_start = next_start+nperseg-noverlap\n",
    "        \n",
    "    # Do row-wise welch\n",
    "    for k in range(col_len):\n",
    "        ff, welch_array[k, :] = signal.welch(sorted_data[k, :], SampleRate,\n",
    "                                             window = 'hann',\n",
    "                                             nperseg = nperseg,\n",
    "                                             nfft = None, detrend = False,\n",
    "                                             return_onesided = True,\n",
    "                                             scaling = 'density', axis = -1,\n",
    "                                             average = 'median')\n",
    "    # Rebin to remove correlations\n",
    "    redff = np.zeros(len(ff[skip:])//rebinSize)\n",
    "    redWelchArr = np.zeros((col_len, len(redff)))\n",
    "    for ii in range(len(redff)):\n",
    "        startInd = skip + rebinSize * ii\n",
    "        endInd = skip + rebinSize * (ii + 1)\n",
    "        redff[ii] = np.mean(ff[startInd:endInd])\n",
    "        for k in range(col_len):\n",
    "            if average=='mean':\n",
    "                redWelchArr[k, ii] = np.mean(welch_array[k, startInd:endInd])\n",
    "            elif average=='median':\n",
    "                redWelchArr[k, ii] = np.median(welch_array[k, startInd:endInd])\n",
    "    \n",
    "    # Reduce PSD to frequencies of interest\n",
    "    redWelchArr = redWelchArr[:, redff>lowerCutoffFreq]\n",
    "    redff = redff[redff>lowerCutoffFreq]\n",
    "    redWelchArr = redWelchArr[:, redff<upperCutoffFreq]\n",
    "    redff = redff[redff<upperCutoffFreq]\n",
    "    \n",
    "    return redff, redWelchArr\n",
    "\n",
    "\n",
    "def logWAtoPSDDist(ff, logWelchArray, nbins=10):\n",
    "    psd_dist = []\n",
    "    for ii, f in enumerate(ff):\n",
    "        logWA = logWelchArray[:, ii]\n",
    "        hist, bin_edges = np.histogram(logWA, bins=nbins)\n",
    "        hist = hist/np.max(hist)\n",
    "        bin_centers = 0.5*(bin_edges[1:] + bin_edges[0:-1])\n",
    "        Q1 = np.percentile(logWA, 25)\n",
    "        Q2 = np.percentile(logWA, 50)\n",
    "        Q3 = np.percentile(logWA, 75)\n",
    "        skewEst = (Q1 + Q2 - 2*Q3)/(Q3 - Q1)\n",
    "        popt, pcov = curve_fit(skewgaus, bin_centers, hist, p0=[np.mean(logWA),np.std(logWA), skewEst, 1],\n",
    "                               maxfev=10000)\n",
    "        psd_dist += [[f, popt, logWA, hist, bin_centers]]\n",
    "    return psd_dist\n",
    "\n",
    "        \n",
    "def tsDatatoPSDDist(tsData, noverlap=None, timeSeg=5, rebinSize=5, skip=3,\n",
    "                    lowerCutoffFreq=70, upperCutoffFreq=600, nbins=10,\n",
    "                    average='median'):\n",
    "    ff, WelchArray = tsDatatoWelchArray(tsData, noverlap=noverlap,\n",
    "                                        timeSeg=timeSeg, rebinSize=rebinSize,\n",
    "                                        skip=skip, lowerCutoffFreq=lowerCutoffFreq,\n",
    "                                        upperCutoffFreq=upperCutoffFreq,\n",
    "                                        average=average)\n",
    "    return logWAtoPSDDist(ff, np.log(WelchArray), nbins)\n",
    "\n",
    "\n",
    "def tsFiletoPSDDist(tsfile, noverlap=None, timeSeg=5, rebinSize=5, skip=3,\n",
    "                    lowerCutoffFreq=70, upperCutoffFreq=600, nbins=10,\n",
    "                    average='median', redo=False):\n",
    "    if isinstance(tsfile, np.ndarray):\n",
    "        tsData = tsfile\n",
    "    elif tsfile.find('TimeSeries') != -1:\n",
    "        tsfile = eu(tsfile)\n",
    "        psdDistFile = tsfile.replace('TimeSeries', 'PSDDist')\n",
    "        if os.path.exists(psdDistFile) and not redo:\n",
    "            print(\"Found calculated PSD Distribution. Reading that...\")\n",
    "            psdDistData = np.loadtxt(psdDistFile)\n",
    "            psd_dist = []\n",
    "            for ii in range(np.shape(psdDistData)[0]):\n",
    "                psd_dist += [[psdDistData[ii, 0], list(psdDistData[ii, 1:]) ]]\n",
    "            return psd_dist\n",
    "        print('Reading time series data...')\n",
    "        tsData = np.loadtxt(tsfile)\n",
    "        print('Done.')\n",
    "    else:\n",
    "        raise RuntimeError('Not a time series file')\n",
    "    print('Fitting PSD Data Distribution...')\n",
    "    psd_dist_1 = tsDatatoPSDDist(tsData, noverlap=noverlap, timeSeg=timeSeg,\n",
    "                                 lowerCutoffFreq=lowerCutoffFreq, upperCutoffFreq=101,\n",
    "                                 rebinSize=rebinSize, skip=skip, nbins=nbins,\n",
    "                                 average=average)\n",
    "    psd_dist_2 = tsDatatoPSDDist(tsData, noverlap=noverlap, timeSeg=timeSeg/10,\n",
    "                                 lowerCutoffFreq=101, upperCutoffFreq=upperCutoffFreq,\n",
    "                                 rebinSize=rebinSize, skip=skip, nbins=nbins,\n",
    "                                 average=average)\n",
    "    print('Done')\n",
    "    psd_dist = psd_dist_1 + psd_dist_2\n",
    "    if isinstance(tsfile, str):\n",
    "        if not os.path.exists(psdDistFile):\n",
    "            print('Storing PSD Distribution Data for later use...')\n",
    "            psdDistData = np.zeros((len(psd_dist), 5))\n",
    "            for ii in range(len(psd_dist)):\n",
    "                psdDistData[ii, 0] = psd_dist[ii][0]\n",
    "                psdDistData[ii, 1] = psd_dist[ii][1][0]\n",
    "                psdDistData[ii, 2] = psd_dist[ii][1][1]\n",
    "                psdDistData[ii, 3] = psd_dist[ii][1][2]\n",
    "                psdDistData[ii, 4] = psd_dist[ii][1][3]\n",
    "            np.savetxt(psdDistFile, psdDistData)\n",
    "    return psd_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Supporting function for likelihood calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bulk\n",
    "def S_Bk_red(freq, nosbud):\n",
    "    S_Bk_num = (4*scc.Boltzmann*nosbud.temp*nosbud.coatStack.WaveLength\n",
    "                        * (1 - nosbud.coatStack.Poisson\n",
    "                           - 2*nosbud.coatStack.Poisson**2))\n",
    "    S_Bk_den = (3*np.pi*freq*nosbud.coatStack.Young\n",
    "                        * ((1 - nosbud.coatStack.Poisson)**2)\n",
    "                        * nosbud.coatStack.Aeff)\n",
    "            # Leaving out substrate from return value\n",
    "    return (S_Bk_num/S_Bk_den)[0:-1]\n",
    "\n",
    "# Shear\n",
    "def S_Sk_red(freq, nosbud):\n",
    "    S_Sk_num = (4*scc.Boltzmann*nosbud.temp*nosbud.coatStack.WaveLength\n",
    "                        * (1 - nosbud.coatStack.Poisson\n",
    "                           - 2*nosbud.coatStack.Poisson**2))\n",
    "    S_Sk_den = (3*np.pi*freq*nosbud.coatStack.Young\n",
    "                        * ((1 - nosbud.coatStack.Poisson)**2)\n",
    "                        * nosbud.coatStack.Aeff)\n",
    "            # Leaving out substrate from return value\n",
    "    return (S_Sk_num/S_Sk_den)[0:-1]\n",
    "\n",
    "def BulkContSlope(freq, nosbud):\n",
    "    if isinstance(freq, Iterable):\n",
    "        return unp.nominal_values(np.array([BulkContSlope(f, nosbud) for f in freq]))\n",
    "    return nosbud.nom*(nosbud.fConv**2)*np.sum(nosbud.q_Bk*S_Bk_red(freq, nosbud))\n",
    "\n",
    "def ShearContSlope(freq, nosbud):\n",
    "    if isinstance(freq, Iterable):\n",
    "        return unp.nominal_values(np.array([ShearContSlope(f, nosbud) for f in freq]))\n",
    "    return nosbud.nom*(nosbud.fConv**2)*np.sum(nosbud.q_Sk*S_Sk_red(freq, nosbud))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Srest(freq, nosbud):\n",
    "    rest = ['coatTO', 'subBr', 'subTE', 'pdhShot',\n",
    "            'pllOsc', 'pllReadout', 'seismic', 'photoThermal','resNPRO']\n",
    "    S_rest = np.zeros_like(freq)\n",
    "    for psd in rest:\n",
    "        PSDest = unp.nominal_values(nosbud.PSDList[psd][0])\n",
    "        S_rest = S_rest + np.interp(freq, nosbud.freq, PSDest)\n",
    "    return S_rest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rem60HzHarm(psd_dist, remNeighbors=False):\n",
    "    ff = np.array([ele[0] for ele in psd_dist])\n",
    "    remInd = []\n",
    "    # Remove any repeated index\n",
    "    for ii in range(1, len(ff)):\n",
    "        if ff[ii] == ff[ii-1]:\n",
    "            remInd = list(set(remInd + [ii]))\n",
    "            \n",
    "    # Remove bad region known to have peaks due to RIN\n",
    "    for ii in range(1, len(ff)):\n",
    "        if ff[ii] > 260 and ff[ii] < 290:\n",
    "            remInd = list(set(remInd + [ii]))\n",
    "    \n",
    "    # Remove 60 Hz harmonics and neighbours\n",
    "    sixtyHarm = np.arange(60, 1000, 60)\n",
    "    for har in sixtyHarm:\n",
    "        if har > ff.min() and har < ff.max():\n",
    "            if np.min(np.abs(ff - har)) < 1:\n",
    "                closestInd = np.argmin(np.abs(ff - har))\n",
    "                remInd = list(set(remInd + [closestInd]))\n",
    "                if closestInd < len(ff)-1:\n",
    "                    remInd = list(set(remInd + [closestInd + 1]))\n",
    "                if closestInd > 0:\n",
    "                    remInd = list(set(remInd + [closestInd - 1]))\n",
    "    \n",
    "    for ind in sorted(remInd, reverse=True):\n",
    "        del psd_dist[ind]\n",
    "    return psd_dist\n",
    "\n",
    "\n",
    "def skewgaus(x, x0, sigma, skewness, prefac):\n",
    "    return prefac*skewnorm.pdf(x, skewness, loc=x0, scale=sigma)\n",
    "\n",
    "\n",
    "def loggaus(x, mean, std):\n",
    "    return -0.5 * ((x - mean) / std)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log-likelihoog calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loglikelihood(BulkLA, ShearLA, BNfile, nosbud, useBNfileOnly=False, overallCompPerc=None):\n",
    "    X, Y = np.meshgrid(BulkLA, ShearLA)\n",
    "    loglkhDist = np.zeros_like(X)\n",
    "    tsfile = BNfile.replace('Spectrum', 'TimeSeries')\n",
    "    psdDistFile = BNfile.replace('Spectrum', 'PSDDist')\n",
    "    if ((os.path.exists(tsfile) or os.path.exists(psdDistFile)) and not useBNfileOnly):\n",
    "        psd_dist = tsFiletoPSDDist(tsfile)\n",
    "        psd_dist = rem60HzHarm(psd_dist)\n",
    "        ff = np.array([ele[0] for ele in psd_dist])\n",
    "        S_rest = Srest(ff, nosbud)\n",
    "        X_B = BulkContSlope(ff, nosbud)\n",
    "        X_S = ShearContSlope(ff, nosbud)\n",
    "        ct = 0\n",
    "        perc = 0\n",
    "        for ii in range(len(ShearLA)):\n",
    "            for jj in range(len(BulkLA)):\n",
    "                estArr = np.log(S_rest + X_B * BulkLA[jj] + X_S * ShearLA[ii])\n",
    "                tempLogLkh = np.zeros_like(ff)\n",
    "                for kk, element in enumerate(psd_dist):\n",
    "                    tempLogLkh[kk] += np.log(skewgaus(estArr[kk], *element[1]))\n",
    "                infPresent = (np.sum(tempLogLkh) == -np.inf)\n",
    "                while infPresent:\n",
    "                    for kk in range(len(ff)):\n",
    "                        if tempLogLkh[kk] == -np.inf:\n",
    "                            if kk > 0 and kk < len(ff) - 1:\n",
    "                                tempLogLkh[kk] = 0.5 * (tempLogLkh[kk-1] + tempLogLkh[kk+1])\n",
    "                            elif kk == 0:\n",
    "                                tempLogLkh[kk] = tempLogLkh[kk+1]\n",
    "                            elif kk == len(ff) - 1:\n",
    "                                tempLogLkh[kk] = tempLogLkh[kk-1]\n",
    "                    infPresent = (np.sum(tempLogLkh) == -np.inf)\n",
    "                loglkhDist[ii, jj] = np.sum(tempLogLkh)\n",
    "                ct = ct + 1\n",
    "                lastperc = perc\n",
    "                perc = np.round(ct*100/len(ShearLA)/len(BulkLA))\n",
    "                if perc != lastperc:\n",
    "                    clear_output()\n",
    "                    if overallCompPerc is not None:\n",
    "                        print('Overall Progress: {}% Completed'.format(str(np.round(overallCompPerc))))\n",
    "                    print('This calculation: {}% Completed'.format(perc))\n",
    "                    sys.stdout.flush()\n",
    "        return loglkhDist\n",
    "    else:\n",
    "        print('Using Beatnote Spectrum file...')\n",
    "        data = np.loadtxt(BNfile)\n",
    "        ff = data[:, 0]\n",
    "        measLogPSD = np.log(data[:, 1]**2)\n",
    "        measLogPSDstd = (np.log(data[:, 3]**2) - np.log(data[:, 2]**2))/2\n",
    "        ff, measLogPSD, measLogPSDstd = cleanData(ff, measLogPSD, measLogPSDstd)\n",
    "        S_rest = Srest(ff, nosbud)\n",
    "        X_B = BulkContSlope(ff, nosbud)\n",
    "        X_S = ShearContSlope(ff, nosbud)\n",
    "        ct = 0\n",
    "        perc = 0\n",
    "        for ii in range(len(ShearLA)):\n",
    "            for jj in range(len(BulkLA)):\n",
    "                estArr = np.log(S_rest + X_B * BulkLA[jj] + X_S * ShearLA[ii])\n",
    "                loglkhDist[ii, jj] = np.sum(loggaus(estArr, measLogPSD, measLogPSDstd))\n",
    "                ct = ct + 1\n",
    "                lastperc = perc\n",
    "                perc = np.round(ct*100/len(ShearLA)/len(BulkLA))\n",
    "                if perc != lastperc:\n",
    "                    clear_output()\n",
    "                    if overallCompPerc is not None:\n",
    "                        print('Overall Progress: {}% Completed'.format(str(np.round(overallCompPerc))))\n",
    "                    print('This calculation: {}% Completed'.format(perc))\n",
    "                    sys.stdout.flush()\n",
    "        return loglkhDist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sweep through the measured data\n",
    "\n",
    "Calculate integrated noise for the data from median value beatnote spectrum calculated during measurement.\n",
    "Take top 10 lowest noise measurements and run ful bayesian analysis on them and choose the measurement that gives lowest most likely bulk loss angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cleanData(ff, beat, beatstd, lowerFreqCutOff=70, upperFreqCutOff=600):\n",
    "    '''\n",
    "    Remove 60 Hz Harmonics and neighbouring bins\n",
    "    '''\n",
    "    beat = beat[ff < upperFreqCutOff]\n",
    "    beatstd = beatstd[ff < upperFreqCutOff]\n",
    "    ff = ff[ff < upperFreqCutOff]\n",
    "    beat = beat[ff > lowerFreqCutOff]\n",
    "    beatstd = beatstd[ff > lowerFreqCutOff]\n",
    "    ff = ff[ff > lowerFreqCutOff]\n",
    "    remInd = []\n",
    "    # Remove any repeated index\n",
    "    for ii in range(1, len(ff)):\n",
    "        if ff[ii] == ff[ii-1]:\n",
    "            remInd = list(set(remInd + [ii]))\n",
    "            \n",
    "    # Remove bad region\n",
    "    for ii in range(1, len(ff)):\n",
    "        if ff[ii] > 260 and ff[ii] < 290:\n",
    "            remInd = list(set(remInd + [ii]))\n",
    "    \n",
    "    # Remove 60 Hz harmonics and neighbours\n",
    "    sixtyHarm = np.arange(60, 1000, 60)\n",
    "    for har in sixtyHarm:\n",
    "        if har > ff.min() and har < ff.max():\n",
    "            closestInd = np.argmin(np.abs(ff - har))\n",
    "            remInd = list(set(remInd + [closestInd]))\n",
    "            if closestInd < len(ff)-1:\n",
    "                remInd = list(set(remInd + [closestInd + 1]))\n",
    "            if closestInd > 0:\n",
    "                remInd = list(set(remInd + [closestInd - 1]))\n",
    "    return np.delete(ff, remInd), np.delete(beat, remInd), np.delete(beatstd, remInd)\n",
    "\n",
    "def updateFromTransRIN(nosbud, BNfile):\n",
    "    transRINfile = BNfile.replace('Spectrum', 'TransRIN')\n",
    "    if os.path.exists(transRINfile):\n",
    "        RINdata = np.loadtxt(transRINfile)\n",
    "        with open(transRINfile, 'r') as f:\n",
    "            header = f.readline()\n",
    "        temp = header[header.find('North DC Val:')+13:].replace(' ', '')\n",
    "        NDC = float(temp.split('Volts')[0])\n",
    "        temp = header[header.find('South DC Val:')+13:].replace(' ', '')\n",
    "        SDC = float(temp.split('Volts')[0])\n",
    "        cpb = int((np.shape(RINdata)[1]-1)/2)\n",
    "        RINdata[:, 1:cpb+1] = RINdata[:, 1:cpb+1]/NDC\n",
    "        RINdata[:, cpb+1:] = RINdata[:, cpb+1:]/SDC\n",
    "\n",
    "        NFin = uf(16700, 1400)\n",
    "        SFin = uf(15100, 340)\n",
    "        PincN = 3.9 * NDC * 1e-3    # Using calibration from 3/20/20\n",
    "        PincS = 4.005 * SDC * 1e-3  # Using calibration from 3/16/20\n",
    "        PcircN = PincN * NFin / np.pi\n",
    "        PcircS = PincS * SFin / np.pi\n",
    "        print('North: Incident Power: {:.2e} mW, Circulating Power: {:.2f} W'.format(PincN*1e3, PcircN))\n",
    "        print('South: Incident Power: {:.2e} mW, Circulating Power: {:.2f} W'.format(PincS*1e3, PcircS))\n",
    "        \n",
    "        nosbud.updatePhotoThermalNoise(RINdata, uf(6, 1) * 1e-6, Pcirc=[PcircN, PcircS])\n",
    "        nosbud.updatePDHShotNoise([PincN, PincS])\n",
    "    return nosbud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# allMeasurements = {}\n",
    "# #''' Uncomment if you have all files and wish to scan them all.\n",
    "# dataDir = ['~/Git/cit_ctnlab/ctn_labdata/data/20200313_SuperBNMeasurement/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200420_SuperBNMeasurement/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200511_SuperBNMeasurement/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200521_SuperBNMeasurement/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200523_SuperBNMeasurementSR/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200528_SuperBNMeasurementSR/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200602_SuperBNMeasurementSR/Data/',\n",
    "#            '~/Git/cit_ctnlab/ctn_labdata/data/20200610_SuperBNMeasurementSR/']\n",
    "#            #'~/Git/cit_ctnlab/ctn_noisebudget/ScienceRun/20200522/',\n",
    "#            #'~/Git/cit_ctnlab/ctn_noisebudget/ScienceRun/20200527/',\n",
    "#            #'~/Git/cit_ctnlab/ctn_noisebudget/Data/dailyBeatNoteData/']\n",
    "# subDir = []\n",
    "# for dr in dataDir:\n",
    "#     for fn in os.listdir(eu(dr)):\n",
    "#         subdr = os.path.join(eu(dr), fn)\n",
    "#         if os.path.isdir(subdr) and fn.find('Data')==0:\n",
    "#             subDir += [subdr]\n",
    "# dataDir += subDir\n",
    "# dataDir = [eu(dr) for dr in dataDir]\n",
    "# for direc in dataDir:\n",
    "#     fl = [fn for fn in os.listdir(direc) if fn.find('Spectrum')!=-1]\n",
    "#     for fn in fl:\n",
    "#         data = np.loadtxt(os.path.join(direc, fn))\n",
    "#         ffclean, beatclean, xx = cleanData(data[:, 0], data[:, 1], data[:, 1])\n",
    "#         intNoise = np.sum(beatclean)\n",
    "#         if intNoise == 0:\n",
    "#             intNoise = np.inf\n",
    "#         allMeasurements[os.path.join(direc, fn)] = intNoise\n",
    "# #'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# def takeEle(x):\n",
    "#     return allMeasurements[x]\n",
    "# lowestNoiseFiles = list(allMeasurements.keys())\n",
    "# lowestNoiseFiles.sort(key=takeEle)\n",
    "# lowestNoiseFiles = lowestNoiseFiles[:60]\n",
    "# for ii, fn in enumerate(lowestNoiseFiles):\n",
    "#     fn = '~' + fn[fn.find('/Git'):]\n",
    "#     lowestNoiseFiles[ii] = fn\n",
    "#     print('----------------------------------------------------------------------------')\n",
    "#     print('Rank:', ii+1)\n",
    "#     print(os.path.basename(fn))\n",
    "#     tsFiletoPSDDist(fn.replace('Spectrum', 'TimeSeries'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if os.uname()[1] == 'ws1':\n",
    "#     from send2trash import send2trash\n",
    "#     todayDir = time.strftime('/home/controls/Git/cit_ctnlab/ctn_labdata/data/20200610_SuperBNMeasurementSR/Data%m%d')\n",
    "#     for dirToClean in dataDir:\n",
    "#         if dirToClean != todayDir:\n",
    "#             allFiles = [os.path.join(dirToClean, fn) for fn in os.listdir(dirToClean) if fn.find('TimeSeries') != -1]\n",
    "#             for fn in allFiles:\n",
    "#                 print(fn)\n",
    "#                 send2trash(fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lowNoiseMeasLog = 'lowNoiseMeas.pkl'\n",
    "SavedPSDsDir = eu('~/Git/cit_ctnlab/ctn_noisebudget/Data/SavedPSDs_20220823/')\n",
    "SavedPSDs = eu('~/Git/cit_ctnlab/ctn_noisebudget/Data/SavedPSDs_20220823/SavedPSDs_20220823_213958.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cwd = eu('~/Git/cit_ctnlab/ctn_noisebudget/BayesianAnalysis')\n",
    "os.chdir(SavedPSDsDir)\n",
    "nosbud = noiseBudget(params='CTN_Noise_Budget_Diff_Loss_Angles.yml')\n",
    "os.chdir(cwd)\n",
    "nosbud.calculateCoatingBrownianNoise();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def updateLog(lowNoiseMeasLog=lowNoiseMeasLog, nosbud=nosbud, SavedPSDs=SavedPSDs):\n",
    "    '''\n",
    "    Redo likelihood calculation for all pickled measurements.\n",
    "    '''\n",
    "    with open(lowNoiseMeasLog, 'rb') as p:\n",
    "        lowNoiseMeas = pickle.load(p)\n",
    "    \n",
    "    minmlBLA = np.inf\n",
    "    BulkLA = lowNoiseMeas['BulkLA']\n",
    "    ShearLA = lowNoiseMeas['ShearLA']\n",
    "    priorDist = logprior(BulkLA, ShearLA)\n",
    "    \n",
    "    for fnind, fn in enumerate(lowNoiseMeas):\n",
    "        overallCompPerc = 100 * fnind / len(lowNoiseMeas)\n",
    "        if fn.find('Spectrum') != -1:\n",
    "            nosbud.loadPSD(SavedPSDs,\n",
    "                           overridePresentFreq=True,\n",
    "                           overridePresentPSD=True)\n",
    "            nosbud = updateFromTransRIN(nosbud, eu(fn))\n",
    "            nosbud.calculateTotalEstNoise()\n",
    "            lkhDist = loglikelihood(BulkLA, ShearLA, eu(fn), nosbud, overallCompPerc=overallCompPerc)\n",
    "            bayProbDist = priorDist + lkhDist\n",
    "            mlBLAind = np.argmax(bayProbDist[0, :])\n",
    "            mlBLA = BulkLA[mlBLAind]\n",
    "            totArea = np.sum(np.exp(bayProbDist[0, :]))\n",
    "            for ii in range(1, len(BulkLA)//2):\n",
    "                if np.sum(np.exp(bayProbDist[0, mlBLAind-ii:mlBLAind+ii])) > 0.9*totArea:\n",
    "                    BLAci90ll = BulkLA[mlBLAind - ii + 1]\n",
    "                    BLAci90ul = BulkLA[mlBLAind + ii - 1]\n",
    "                    BLAci90pm = (BLAci90ul - BLAci90ll)/2\n",
    "                    break\n",
    "            lowNoiseMeas[fn] = {'mlBLA': mlBLA,\n",
    "                                'BLAci90pm': BLAci90pm,\n",
    "                                'lkhDist': lkhDist}\n",
    "            if mlBLA < minmlBLA:\n",
    "                BNfile = fn\n",
    "                minmlBLA = mlBLA\n",
    "                lowNoiseMeas['minmlBLA'] = minmlBLA\n",
    "                lowNoiseMeas['BNfile'] = BNfile\n",
    "    with open(lowNoiseMeasLog, 'wb') as logFile:\n",
    "        pickle.dump(lowNoiseMeas, logFile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# updateLog()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    with open(lowNoiseMeasLog, 'rb') as p:\n",
    "        lowNoiseMeas = pickle.load(p)\n",
    "except BaseException:\n",
    "    BulkLA = np.arange(0, 16e-4, 1e-6)\n",
    "    ShearLA = np.array([5.2e-7])\n",
    "    priorDist = logprior(BulkLA, ShearLA)\n",
    "    lowNoiseMeas = {'BulkLA': BulkLA,\n",
    "                    'ShearLA': ShearLA,\n",
    "                    'priorDist': priorDist}\n",
    "try:\n",
    "    minmlBLA = lowNoiseMeas['minmlBLA']\n",
    "except BaseException:\n",
    "    minmlBLA = np.inf\n",
    "        \n",
    "BulkLA = lowNoiseMeas['BulkLA']\n",
    "ShearLA = lowNoiseMeas['ShearLA']\n",
    "priorDist = logprior(BulkLA, ShearLA)\n",
    "\n",
    "try:\n",
    "    for fnind, fn in enumerate(lowestNoiseFiles):\n",
    "        overallCompPerc = 100 * fnind / len(lowNoiseMeas)\n",
    "        if fn not in lowNoiseMeas:\n",
    "            nosbud.loadPSD(SavedPSDs,\n",
    "                           overridePresentFreq=True, overridePresentPSD=True)\n",
    "            nosbud = updateFromTransRIN(nosbud, eu(fn))\n",
    "            nosbud.calculateTotalEstNoise()\n",
    "            lkhDist = loglikelihood(BulkLA, ShearLA, eu(fn), nosbud, overallCompPerc=overallCompPerc)\n",
    "            bayProbDist = priorDist + lkhDist\n",
    "            mlBLAind = np.argmax(bayProbDist[0, :])\n",
    "            mlBLA = BulkLA[mlBLAind]\n",
    "            totArea = np.sum(np.exp(bayProbDist[0, :]))\n",
    "            for ii in range(1, len(BulkLA)//2):\n",
    "                if np.sum(np.exp(bayProbDist[0, mlBLAind-ii:mlBLAind+ii])) > 0.9*totArea:\n",
    "                    BLAci90ll = BulkLA[mlBLAind - ii + 1]\n",
    "                    BLAci90ul = BulkLA[mlBLAind + ii - 1]\n",
    "                    BLAci90pm = (BLAci90ul - BLAci90ll)/2\n",
    "                    break\n",
    "            lowNoiseMeas[fn] = {'mlBLA': mlBLA,\n",
    "                                'BLAci90pm': BLAci90pm,\n",
    "                                'lkhDist': lkhDist}\n",
    "            if mlBLA < minmlBLA:\n",
    "                BNfile = fn\n",
    "                minmlBLA = mlBLA\n",
    "                lowNoiseMeas['minmlBLA'] = minmlBLA\n",
    "                lowNoiseMeas['BNfile'] = BNfile\n",
    "            with open(lowNoiseMeasLog, 'wb') as logFile:\n",
    "                pickle.dump(lowNoiseMeas, logFile)\n",
    "except BaseException:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BNfile = lowNoiseMeas['BNfile']\n",
    "print('Best measurement found:')\n",
    "print(os.path.basename(BNfile))\n",
    "mlBLA = lowNoiseMeas[BNfile]['mlBLA']\n",
    "BLAci90pm = lowNoiseMeas[BNfile]['BLAci90pm']\n",
    "print('Estimated Bulk Loss Angle {:.2e} +- {:.2e} radians assuming logPSD is normal'.format(mlBLA, BLAci90pm))\n",
    "lkhDist = lowNoiseMeas[BNfile]['lkhDist']\n",
    "bayProbDist = priorDist + lkhDist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BLAstr = r'\\Phi_\\text{B} = ' + '({:.2f} \\pm {:.2f})'.format(mlBLA*1e4, BLAci90pm*1e4) + r' \\times 10^{-4}'\n",
    "print(BLAstr)\n",
    "#with open(eu('~/Git/cit_ctnlab/ctn_paper/data/Bulk_Loss_Fit_Value_String.tex'), 'w') as f:\n",
    "#    f.writelines(BLAstr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bayesFactor(model1='lowNoiseMeas.pkl', model2='lowNoiseMeasWithSlope.pkl'):\n",
    "    with open(model1, 'rb') as p:\n",
    "        model1 = pickle.load(p)\n",
    "    minmlBLA = np.inf\n",
    "    for key, value in model1.items():\n",
    "        if key.find('.txt')!=-1:\n",
    "            if value['mlBLA'] < minmlBLA:\n",
    "                minmlBLA = value['mlBLA']\n",
    "                minKey = key\n",
    "    dBLA = model1['BulkLA'][1] - model1['BulkLA'][0]\n",
    "    evidence1 = np.sum(np.exp(model1['priorDist'] + model1[minKey]['lkhDist']) * dBLA)\n",
    "    print('Evidence for model 1 is {:.2e}'.format(evidence1))\n",
    "    \n",
    "    with open(model2, 'rb') as p:\n",
    "        model2 = pickle.load(p)\n",
    "    minmlBLA = np.inf\n",
    "    for key, value in model2.items():\n",
    "        if key.find('.txt')!=-1:\n",
    "            if value['mlBLA'] < minmlBLA:\n",
    "                minmlBLA = value['mlBLA']\n",
    "                minKey = key\n",
    "    dBLA = model2['BulkLA'][1] - model2['BulkLA'][0]\n",
    "    dBLAslope = model2['BulkLAslope'][1] - model2['BulkLAslope'][0]\n",
    "    evidence2 = np.sum(np.exp(model2['priorDist'] + model2[minKey]['lkhDist']) * dBLA * dBLAslope)\n",
    "    print('Evidence for model 2 is {:.2e}'.format(evidence2))\n",
    "    \n",
    "    print('Bayes factor for model2:mode1 is {:.2e}'.format(evidence2/evidence1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayesFactor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[16,12])\n",
    "ax = fig.gca()\n",
    "exppriorDist = np.exp(priorDist[0, :])/np.max(np.exp(priorDist[0, :]))\n",
    "\n",
    "color = 'tab:blue'\n",
    "ax.plot(BulkLA * 1e4, exppriorDist, color=color,\n",
    "        lw=4, label='Likelihood Probability (Norm.)')\n",
    "ax.set_xlabel('Bulk Loss Angle (' + r'$\\times 10^{-4}$' + ' rad)')\n",
    "ax.set_ylabel('Prior Probability (normalized)',\n",
    "              color=color)\n",
    "ax.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "ax2 = ax.twinx()\n",
    "color='tab:orange'\n",
    "ax2.plot(BulkLA * 1e4, priorDist[0, :], color=color,\n",
    "         lw=4, label='Log-Prior Distribution (Unnorm.)')\n",
    "ax2.set_ylabel('Log-Prior Dist. (unnormalized)',\n",
    "               color=color)\n",
    "ax2.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "ax.set_title('Prior Probability Distribution of Bulk loss angles')\n",
    "figlist = [fig]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[16,12])\n",
    "ax = fig.gca()\n",
    "explkhDist = np.exp(lkhDist[0, :])/np.max(np.exp(lkhDist[0, :]))\n",
    "\n",
    "\n",
    "color = 'tab:blue'\n",
    "ax.plot(BulkLA * 1e4, explkhDist, color=color,\n",
    "        lw=4, label='Likelihood Probability (Norm.)')\n",
    "ax.set_xlabel('Bulk Loss Angle (' + r'$\\times 10^{-4}$' + ' rad)')\n",
    "ax.set_ylabel('Likelihood Probability (normalized)',\n",
    "              color=color)\n",
    "ax.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "ax2 = ax.twinx()\n",
    "color='tab:orange'\n",
    "ax2.plot(BulkLA * 1e4, lkhDist[0, :], color=color,\n",
    "         lw=4, label='Log-Likelihood Distribution (Unnorm.)')\n",
    "ax2.set_ylabel('Log-Likelihood Dist. (unnormalized)',\n",
    "               color=color)\n",
    "ax2.tick_params(axis='y', labelcolor=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig.savefig(eu('~/Git/cit_ctnlab/ctn_paper/figures/LikelihoodProbDist.pdf'),\n",
    "#            facecolor=fig.get_facecolor(),\n",
    "#            bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig.savefig(eu('~/Git/writing-presenting/candidacyTalk/figures/LikelihoodProbDist.pdf'),\n",
    "#            facecolor=fig.get_facecolor(),\n",
    "#            bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.set_title('Likelihood Distribution of Bulk loss angles for'\n",
    "             'measured ASD of beatnote\\n from 70 Hz to 600 Hz')\n",
    "figlist += [fig]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[16,12])\n",
    "ax = fig.gca()\n",
    "expbayProbDist = np.exp(bayProbDist[0, :])/np.max(np.exp(bayProbDist[0, :]))\n",
    "\n",
    "\n",
    "color = 'tab:blue'\n",
    "ax.plot(BulkLA * 1e4, expbayProbDist, color=color,\n",
    "        lw=4, label='Bayesian Probability (Norm.)')\n",
    "ax.set_xlabel('Bulk Loss Angle (' + r'$\\times 10^{-4}$' + ' rad)')\n",
    "ax.set_ylabel('Probability Distribution (normalized)',\n",
    "              color=color)\n",
    "ax.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "ax2 = ax.twinx()\n",
    "color='tab:orange'\n",
    "ax2.plot(BulkLA * 1e4, bayProbDist[0, :], color=color,\n",
    "         lw=4, label='Log-Bayesian Probability Distribution (Unnorm.)')\n",
    "ax2.set_ylabel('Log-Prob. Dist. (unnormalized)',\n",
    "               color=color)\n",
    "ax2.tick_params(axis='y', labelcolor=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.set_title('Bayesian Inferred Porbability Distribution of Bulk loss angles for'\n",
    "             'measured ASD of beatnote\\nfrom 70 Hz to 600 Hz')\n",
    "figlist += [fig]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayProbData = np.zeros((len(BulkLA), 4))\n",
    "bayProbData[:, 0] = BulkLA\n",
    "bayProbData[:, 1] = exppriorDist\n",
    "bayProbData[:, 2] = explkhDist\n",
    "bayProbData[:, 3] = expbayProbDist\n",
    "# np.savetxt(eu('~/Git/cit_ctnlab/ctn_paper/'\n",
    "#              '/data/bayProbData.txt'), bayProbData,\n",
    "#           header='Bulk Loss Angle    Prior    Likelihood    Bayesian')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nosbud.loadPSD(SavedPSDs,\n",
    "               overridePresentFreq=True, overridePresentPSD=True)\n",
    "nosbud.loadASD(eu(BNfile),\n",
    "               label='Measured Beatnote Spectrum',\n",
    "               key='beat')\n",
    "nosbud = updateFromTransRIN(nosbud, eu(BNfile))\n",
    "nosbud.calculateTotalEstNoise()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,