# Here we make the list of the free models
# the fixed component is handled seperately
modelList = [ draco, galdif ] 
par_index = 0 # This is the index of the source we care about (i.e., Draco)

"""
Explanation: Ok, so we have 23616 counts in our ROI, and the best-fit sum of the Galactic diffuse model and all the fixed model components adds up to 23705 counts.
The counts for the Draco model are for a powerlaw spectrum with and index of -2 and an prefactor of $10^{-12} cm^{-2}s^{-1}MeV{-1}$ at 1000 GeV. 
We are going to be fitting for the scale factors of the Draco and Galactic diffuse components with respect to those values.
End of explanation
"""


help(lfu.NLL_func)

ftomin = lfu.NLL_func(n_obs,fixed,modelList)
init_pars = np.ones((2))

print "NLL = %.1f"%ftomin(init_pars)

"""
Explanation: Ok, now we are going to construct a function that will return the negative log-likelihood for a particular set of normailzation of the draco and Galactic Diffuse flux.
There is utility function to do this in the LikeFitUtils.py file in this directory called NLL_func.   In terms of the objects we are using the formula for the negative log-likelihood is:
$-\log\mathcal{L} = \sum_{i} n_{i} \log ( m_{i,\rm{fixed}} + \sum_{j} p_{j} m_{ij} ) - m_{i,\rm{fixed}} - \sum_{j} p_{j} m_{ij}$
Where the index $i$ runs over all the energy bins and pixels, and $j$ run over the model compoents 
* $n_{i}$ are the observed counts: (n_obs)
* $m_{i,\rm{fixed}}$ are the preditcted counts of the fixed model component: (fixed)
* $m_{ij}$ are the predicted counts of the free model components: [draco, galdif]
* $p_{j}$ are the fit parameters: the will be the input to function
End of explanation
"""


help(lfu.Minimize)

# Ok, let's fit the function and parse out the results
result = lfu.Minimize(ftomin,init_pars)
mle_pars = result[0]
mle = mle_pars[par_index]
nll_min = result[1]
nll_null = ftomin([0.,mle_pars[1]])
TS = 2*(nll_null-nll_min)
pvalue = 0.5 # fixme
print "Maximum likelihood estimate = %.1e cm-2 s-1 MeV-1"%(mle*1e-12)
print "Minimum function value = %.2f"%nll_min
print "Test Statistic = %.2f"%TS
print "p-value = %.2f"%pvalue

"""
Explanation: Ok, now we are going to want to minize the function.   I added a function to LikeFitUtils.py to do this using scipy.optimize.fmin, the function is called Minimize
End of explanation
"""


# Set up a likelihood scan
par_bounds = (1e-2,2.0)
nsteps = 25

# Make a set of input parameter vectors that will serve to do the scan
par_sets = lfu.MakeParSets_1DScan(mle_pars,par_index,par_bounds[0],par_bounds[1],nsteps,log=True)
print par_sets

"""
Explanation: Ok, so here we have seen that the best-fit value is $1.0^{-13}$cm$^{-2}$s$^{-1}$MeV$^{-1}$, but that the test statistic is only 0.47.  This tells us that we have not significantly detected emission from Draco with a powerlaw index of -2.
Question:
Why does $TS = 0.47$ imply a non-detection?
In the absence of a detection, we are now going to derive upper limits on the powerlaw prefactor.   To do this are going to have to calculate the likelihood as a function of the Draco flux normalization.
End of explanation
"""


# This tells us to fix the Draco normalization during the Profile likelihood scan
fix_par_mask = np.zeros((2),'?')
fix_par_mask[par_index] = True

pf1 = lfu.ParameterScan(ftomin,par_sets)
print pf1

# Let's plot the likelihood, actually what we are going to plot is the 
# delta log-likelihood (w.r.t. the maximum)
fig1,ax1 = lfu.PlotNLLScan(par_sets[0:,par_index],nll_min-pf1)

"""
Explanation: Now, we are going to scan likelihood two different ways.
1. Keeping the normalization of the other parameters fixed
1. Re-optimizing the normalization of the other paramters are each scan step, this is called the profile likelihood.
End of explanation
"""


interp1 = lfu.BuildInterpolator(par_sets,par_index,pf1)

"""
Explanation: In our case it doesn't really matter because the fit is simple, but it can be useful to interpolate the likelihood between the scan points rather that recompute it.
End of explanation
"""


help(lfu.SolveForErrorLevel)

lim1 = lfu.SolveForErrorLevel(interp1,nll_min,1.35,mle,par_bounds)
print "Simple upper limit %.2e cm-2 s-1 MeV-1"%(lim1*1e-12)

"""
Explanation: Now we are going to solve for the 95% confidence level upper limits.  In the likelihood ratio test where the test hypothesis has one extra degree of freedom, and we are at the boundary of that additional degree of freedom this occurs at the point that the log-likelihood is 1.35 less than maximum. 
We can derive that factor of 1.35 as follows.

Chernoff's theorem states that if the null-hypothesis is correct then 50% of the trials will have TS=0 and the other 50% will be $\chi^2$-distributed with one degree of freedom.
We assign the maximum likelihood estimate as the null-hypothesis.
For 1 degree of freedom the cumulative distribution function of the $\chi^2$-distribution is simply $erf(\sqrt{TS/2})$, and the p-value is complement of that, or $erfc(\sqrt{TS/2})$.
Therefore, in our case the p-value as a function x is $p = 0.5 erfc(\sqrt{TS/2})$.
Recall that TS is twice the delta log-likelihood, so the TS/2 is simply $\Delta\mathcal{L}$, and $p=0.5 erfc(\sqrt{\Delta\mathcal{L}})$.
We can solve numerically for $p=0.05$, or any other confidence level for that matter.
For $p=0.05$ we get $\Delta\mathcal{L} = 1.3528$.

Ok, let's solve for the upper limits now.
End of explanation
"""


fig2,ax2 = lfu.PlotNLLScan(par_sets[0:,par_index],nll_min-pf1)

ax2.hlines(-1.35,0,2.0,linestyles=u'dotted')
ax2.vlines(lim1,-10,1,linestyles=u'dashed')
    
leg = ax2.legend()

"""
Explanation: Ok, now let's make a plot of this.
End of explanation
"""


help(lfu.ProfileScan)

pf2 = lfu.ProfileScan(n_obs, par_sets, fix_par_mask, fixed, modelList)
interp2 = lfu.BuildInterpolator(par_sets,par_index,pf2)
lim2 = lfu.SolveForErrorLevel(interp2,nll_min,1.35,mle,par_bounds)
print "Delta Log-Likelihood at 1.0e-12 cm-2 s-1 MeV-1 is %.1f"%(interp2(1.0)-nll_min)
print "Profile upper limit %.2e cm-2 s-1 MeV-1"%(lim2*1e-12)

# Ok, let's plot both sets of limits
fig3,ax3 = lfu.PlotNLLScan(par_sets[0:,par_index],nll_min-pf1)
ax3.plot(par_sets[0:,par_index],nll_min-pf2,'b-',label="Profile")

ax3.hlines(-1.35,0,2.0,linestyles=u'dotted')
ax3.vlines(lim1,-10,1,linestyles=u'dashed')
    
leg = ax3.legend()

"""
Explanation: The dashed vertical line indicates the upper limit we have computed: 4.68e-13 cm-2 s-1 MeV-1
Profile likelihood
The potential flaw in what we have done is that we didn't allow any of the other paramters to vary as we scanned the log-likelihood as a function of the target flux.  To do things correctly we should at least allow the normalization of the diffuse background to vary to compenstate for the variation in the target flux.
We will do this by re-optimizing the normalization of the other paramters are each scan step, this is called the profile likelihood.
Question, how do we expect this to change our results?
End of explanation
"""


