graph.BPT(df)

"""
Explanation: Bayesian Inference of Censored SDSS data

By: Christopher Stewart
Overview


Introduce the problem
Data Selection
Issues with traditional approach
Introduce the BPT
Limit Cases


Handling Limits
Creating the model
Example Case
Results


What next

Intro  -  Data Selection

Primary SQL quarry

Selected from redshift range of 0.07 to 0.1
error above zero for following lines
$\text{N}_{[\text{II}]}$
$\text{O}_{[\text{III}]}$
$\text{He}_{[\text{II}]}$
$\text{H}_\alpha$
$\text{H}_\beta$



Total Objects: 42239
Sigma Filtering
Objects were removed from the survey following the condition:
$0 <\text{flux_err} < 1\times e^{-15}$
Intro  -  Issues with traditional approach
End of explanation
"""


no_limits = df[(oiii_healthy & nii_healthy)]
graph.BPT(no_limits)

"""
Explanation: No limits in any ratio
End of explanation
"""


print('Number of Objects\n-------------\nAll Data: {}\nWell Def: {}'.format(df.shape[0],no_limits.shape[0]))
print('Loss Percentage: {0:.3f}%'.format((1-no_limits.shape[0]/df.shape[0])*100))

"""
Explanation: How many object we lose if we only take well defined objects
End of explanation
"""


print('OIII upLim {}\nOIII lowLim: {}\nNII upLim: {}\nNII lowLim: {}'.format(df[oiii_ratio_upLim & nii_healthy].shape[0],
                                                                            df[oiii_ratio_lowLim & nii_healthy].shape[0],
                                                                            df[nii_ratio_upLim & oiii_healthy].shape[0],
                                                                            df[nii_ratio_lowLim & oiii_healthy].shape[0]))

"""
Explanation: Lose around 9.6% of our objects
Introduce the BPT

Each object on the BPT hase two ratios which are made up of two fluxes with individual sigma values.
For example a single object has the following data:   
$\log_{10}(\frac{\text{O}{\text{III}}}{\text{H}{\beta}}) \rightarrow \left[
\begin{array}{ll}
      \text{O}{\text{flux}} & \text{O}{{\sigma}} \
      \text{H}\beta_{\text{flux}} & \text{H}\beta_{{\sigma}}  \
\end{array} 
\right]$
$\log_{10}(\frac{\text{N}{\text{II}}}{\text{H}{\alpha}}) \rightarrow \left[
\begin{array}{ll}
      \text{N}{\text{flux}} & \text{N}{{\sigma}} \
      \text{H}\alpha_{\text{flux}} & \text{H}\alpha_{{\sigma}}  \
\end{array} 
\right]$
End of explanation
"""


print(d_limits)

"""
Explanation: Limit cases

Limit:
$\text{flux} < \text{flux error}$
Lower Limit:
Let  $f = \frac{A}{B} \text{ where } B < \sigma_B \text{, so  } f = \frac{A}{\sigma_B} \text{ is an lower limit.}$
Upper Limit:
Let  $f = \frac{A}{B} \text{ where } A < \sigma_A \text{, so } f = \frac{\sigma_A}{B} \text{ is an upper limit.}$
1D Limit
Is a limit in only one ratio. For example a lower limit in x would have the following shape,

2D Limit
Is a limit in each ratio.  For example a lower limit in x and y would have the following shape,

Number of Limit cases

1D Limits:
End of explanation
"""


print(dd_limits)

"""
Explanation: 2D Limits:
End of explanation
"""


test_obj = df.iloc[[21226]]
test_obj[['Log(O_III/H_Beta)','Log(O_III/H_Beta)_Sigma','Log(N_II/H_Alpha)','Log(N_II/H_Alpha)_Sigma', 'O_III/H_Beta', 'H_Beta', 'H_Beta_Sigma', 'O_III']]

graph.BPT(test_obj)
test_obj[lim_columns]

"""
Explanation: Handling Limits

Creating the model
Example Case:
Object has a limit in OIII so this is an upper limit in our y
End of explanation
"""


x_test_data = [test_obj.loc[21226,'Log(N_II/H_Alpha)'],test_obj.loc[21226,'Log(N_II/H_Alpha)_Sigma'],
               test_obj.loc[21226,'N_II'],test_obj.loc[21226,'N_II_Sigma'],
               test_obj.loc[21226,'H_Alpha'],test_obj.loc[21226,'H_Alpha_Sigma']
              ]
y_test_data = [test_obj.loc[21226,'Log(O_III/H_Beta)'],test_obj.loc[21226,'Log(O_III/H_Beta)_Sigma'],
               test_obj.loc[21226,'O_III'],test_obj.loc[21226,'O_III_Sigma'],
               test_obj.loc[21226,'H_Beta'],test_obj.loc[21226,'H_Beta_Sigma']
              ]

@sampled
def oiii_limit(xData,yData):
    x,x_sigma,x_a,x_a_sigma,x_b,x_b_sigma = xData
    y,y_sigma,y_a,y_a_sigma,y_b,y_b_sigma = yData
    
    model_x_a = pm.Normal('N_II', mu=x_a, sd = x_a_sigma)
    model_x_b = pm.Normal('H_Alpha', mu=x_b, sd = x_b_sigma)
    model_y_b = pm.Normal('H_Beta', mu=y_b, sd = y_b_sigma)
    #the upper bound of 0.94 comes from solving ratio OIII/HBeta
    model_y_a_limit = pm.Uniform('OIII_Uniform',lower =0 ,upper=0.9417)
    model_y_a_gauss = pm.Normal('OIII_Normal',mu = 0.9417, sd = y_a_sigma)
    y_use_limit = pm.Binomial('OIII_swap', n=1,p=0.5)
    model_y_a = pm.Normal('OIII', mu = tt.switch(y_use_limit,model_y_a_limit,model_y_a_gauss), sd =1)
    
    model_NII_HAlpha = pm.Deterministic('NII/HAlpha', model_x_a/model_x_b)
    model_OIII_HBeta = pm.Deterministic('OIII/HBeta', model_y_a/model_y_b)
    


"""
Explanation: Creating model in pymc3
End of explanation
"""


with oiii_limit(xData = x_test_data,yData =y_test_data):
    base_trace = pm.sample(10000)

"""
Explanation: Run model
End of explanation
"""


pm.traceplot(base_trace, varnames = ['OIII_Normal','OIII_Uniform'])


def plot_jointkde(x_pdf,y_pdf):
    def StarForming(x): return (1.3 + 0.61 / (x - 0.04))
    g = sns.jointplot(np.log10(x_pdf),np.log10(y_pdf),kind = 'kde',size=10)
    xSpan = np.linspace(-0.6, -0.1, num=200)
    g.ax_joint.plot(xSpan, StarForming(xSpan), 'k--',-0.31,-0.37, 'bo')
    g.set_axis_labels('Log(N_II/H_Aplha)','Log(O[III]/H_Beta)')



"""
Explanation: OIII in depth
End of explanation
"""


plot_jointkde(x_pdf= base_trace['NII/HAlpha'],y_pdf =base_trace['OIII/HBeta'])

pm.summary(base_trace)



"""
Explanation: Our 2D KDE
With our original point plotted as a blue dot
End of explanation
"""


def test_percentAGN(x,x_sig,y,y_sig):
    def AGNLine(x): return 1.19 + 0.61 / (x - 0.47)
    def SFLine(x): return 1.3 + 0.61 / (x - 0.04)
    lineCorrectionAGN = 0.2
    lineCorrectionSF = -0.2

    x0_array = np.float64(np.ones((100,)) * -x)
    y0_array = np.float64(np.ones((100,)) * -y)
    xLineSpace = np.linspace((x - x_sig), (x + x_sig), 100)
    yLineSpace = np.linspace((y - y_sig), (y + y_sig),
                             100)[:, None]
    ellipse = ((xLineSpace - x) / x_sig)**2 + ((yLineSpace - y) / y_sig)**2 <= 1
    boolean_AGN_line = np.subtract(
        yLineSpace, AGNLine(xLineSpace)) > 0  # AGN fit line
    boolean_AGN_Correction = (np.subtract(xLineSpace, np.ones(
        (100, 100)) * lineCorrectionAGN) > 0)  # Fixes problem with the fit line
    # combines fit line and correction into single boolean array
    boolean_AGN_Mask = np.where(np.logical_or(
        boolean_AGN_line, boolean_AGN_Correction), True, False)
    
    # plt.imshow(boolean_AGN_Mask, origin = 'lower') #for debugging

    # Starforming Fit Lines
    boolean_SF_line = np.subtract(yLineSpace, SFLine(
        xLineSpace)) < 0  # Starforming fit line
    boolean_SF_Correction = (np.subtract(xLineSpace, np.ones(
        (100, 100)) * lineCorrectionSF) < 0)  # Fixes problem with the fit line
    # combines fit line and correction into single boolean array
    boolean_SF_Mask = np.where(np.logical_and(
        boolean_SF_line, boolean_SF_Correction), True, False)

    # Neither mask indicates areas that are neither AGN or SF
    boolean_Neither_Mask = np.where(np.logical_or(
        boolean_SF_Mask, boolean_AGN_Mask), False, True)
    # plt.imshow(boolean_SF_Mask, origin = 'lower') #for debugging

    # Creat probability distribution function based on array index locations
    xWeight = (np.exp(-(np.add(xLineSpace, x0_array)**2) *
                      (2 * x_sig**2)) * (2 * math.pi * x_sig**2))**2
    yWeight = (np.exp(-(np.add(yLineSpace, y0_array)**2) *
                      (2 * y_sig**2)) * (2 * math.pi * y_sig**2))**2
    weight_Matrix = np.sqrt(xWeight + yWeight)

    # AGN total weight is the sum of weight matrix elements that are true in both the ellipse boolean mask and agn boolean mask
    AGN_Weight = np.sum(
        weight_Matrix[np.logical_and(ellipse, boolean_AGN_Mask)])
    SF_Weight = np.sum(weight_Matrix[np.logical_and(ellipse, boolean_SF_Mask)])
    Neither_Weight = np.sum(
        weight_Matrix[np.logical_and(ellipse, boolean_Neither_Mask)])
    # total weight of whole ellipse
    Total_Weight = np.sum(weight_Matrix[ellipse])

    # Percenbt AGN or starforming
    AGN_Percent = AGN_Weight / Total_Weight
    SF_Percent = SF_Weight / Total_Weight
    Neither_Percent = Neither_Weight / Total_Weight
    AGN_Percent += Neither_Percent
    
    return AGN_Percent

"""
Explanation: Point location Result

Old:
$$\log_{10}(\frac{\text{O}{\text{III}}}{\text{H}{\beta}}) =   -0.37$$  
$$\log_{10}(\frac{\text{N}{\text{II}}}{\text{H}{\alpha}}) = -0.31 $$
New:
$$\log_{10}(\frac{\text{O}{\text{III}}}{\text{H}{\beta}}) = -0.67$$ 
$$\log_{10}(\frac{\text{N}{\text{II}}}{\text{H}{\alpha}}) = -0.31$$
End of explanation
"""


test_obj[['BPT:P(AGN)']]

"""
Explanation: P(AGN) Result

Old:
End of explanation
"""


test_percentAGN(x=-0.31,x_sig=0.03,y=-0.67,y_sig=np.log10(4.26))

"""
Explanation: New:
End of explanation
"""


