Metadata-Version: 2.1
Name: LensCalcPy
Version: 0.0.2
Summary: Calculate Microlensing Observables
Home-page: https://github.com/NolanSmyth/LensCalcPy
Author: Nolan Smyth
Author-email: nolanwsmyth@gmail.com
License: Apache Software License 2.0
Keywords: nbdev jupyter notebook python
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Developers
Classifier: Natural Language :: English
Classifier: Programming Language :: Python :: 3.7
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: License :: OSI Approved :: Apache Software License
Requires-Python: >=3.7
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: matplotlib
Requires-Dist: astropy
Requires-Dist: pathos
Requires-Dist: tqdm
Provides-Extra: dev

LensCalcPy
================

<!-- WARNING: THIS FILE WAS AUTOGENERATED! DO NOT EDIT! -->

Will write a better readme soon

## Install

``` sh
pip install LensCalcPy
```

## How to use

We can calculate the distribution of crossing times for a given PBH
population

``` python
f_pbh = 1 # fraction of dark matter in PBHs

ts = np.logspace(-2, 1, 30)
pbhs = [Pbh(10**(i), f_pbh) for i in np.linspace(-9, -7, 3)]
result = np.zeros((len(pbhs), len(ts)))
for i, pbh in enumerate(pbhs):
    result[i, :] = pbh.compute_differential_rate(ts)
```

    /Users/nolansmyth/opt/anaconda3/lib/python3.9/site-packages/scipy/integrate/quadpack.py:879: IntegrationWarning: The integral is probably divergent, or slowly convergent.
      quad_r = quad(f, low, high, args=args, full_output=self.full_output,

``` python
for i, pbh in enumerate(pbhs):
    plt.loglog(ts, result[i], label=r"$M_{\rm{PBH}} = $" + scientific_format(pbh.mass,0) + "$M_{\odot}$")

plt.xlabel(r"$t_E$ [h]", fontsize=16)
plt.ylabel(r"$d\Gamma/dt$ [events/star/hr/hr]", fontsize=16)
plt.xlim(1e-2, 1e1)
plt.ylim(1e-10, 1e-4)

plt.legend()
plt.show()
```

![](index_files/figure-commonmark/cell-3-output-1.png)

Similarly, we can calculate the distribution of crossing times for an
FFP population

``` python
alpha = 2
fp = Ffp(alpha)
```

``` python
ts = np.logspace(-2, np.log10(1e3), num=30)
diff_rates = fp.compute_differential_rate(ts, finite=False)
```

``` python
#Note: event rate normalized to 1 FFP/ISO per host star
plt.loglog(ts, diff_rates, color="blue")
plt.xlabel(r"$t_E$ [h]", fontsize=16)
plt.ylabel(r"$d\Gamma/dt$ [events/star/hr/hr]", fontsize=16)
plt.xlim([0.009, 1e3])
plt.ylim([1e-32, 1e-13])
plt.show()
```

![](index_files/figure-commonmark/cell-6-output-1.png)

## Statistics

Eventually will be able to perform statistical analysis on the
distributions of crossing times for a given population and compare with
output from [PopSyCLE](https://github.com/jluastro/PopSyCLE/). For now,
here’s a toy example assuming perfect mass reconstruction

``` python
# Generate example observed counts with a bump

min_bin = 1
max_bin = 10
bin_edges = np.linspace(min_bin,max_bin,20)  # Assuming bins from 1 to 10
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2


a = 200 # normalization (Proxy for overall count of FFPs)
index = -2 # power law index (FFP Power Law Index)
bump_position = 6 # position of bump in x (PBH mass)
bump_height = 20 # extra counts at position of bump (PBH counts)

observed_counts = generate_observed_counts_with_bump(bin_centers, a, index, bump_position, bump_height)
optimized_params = get_MLE_params(bin_centers, observed_counts)

expected_counts_opt = generate_observed_counts_with_bump(bin_centers, *optimized_params)

print("Optimized parameters: a = {:.2f}, index = {:.2f}, bump_position = {:.2f}, bump_height = {:.2f}".format(*optimized_params))
print("True parameters: a = {:.2f}, index = {:.2f}, bump_position = {:.2e}, bump_height = {:.2f}".format(a, index, bump_position, bump_height))
```

    Optimized parameters: a = 217.97, index = -2.10, bump_position = 6.16, bump_height = 16.85
    True parameters: a = 200.00, index = -2.00, bump_position = 6.00e+00, bump_height = 20.00

``` python
# Plot the observed counts and the initial/optimized expected counts
plt.figure(figsize=(6, 4))
# plt.subplot(121)
plt.bar(bin_centers, observed_counts, alpha=0.7, label="Observed", width=np.diff(bin_edges))
plt.plot(bin_centers, power_law(bin_centers, a=a, index=index), label="Null Hypothesis", color='r', linewidth=2)
plt.plot(bin_centers, expected_counts_opt, label="MLE Params", color='g', linewidth=2)

plt.xlabel("Bin center")
plt.ylabel("Counts")
# plt.yscale("log")
plt.legend()
plt.title("Observed Counts vs Initial Expected Counts")


plt.tight_layout()
plt.show()
```

![](index_files/figure-commonmark/cell-8-output-1.png)

``` python
# Calculate the likelihood ratio test statistic
lr_statistic = llr_test(bin_centers, observed_counts)

# Calculate the p-value using a chi-square distribution
p_value = chi2.sf(lr_statistic, df=2)  # Two degrees of freedom for the difference in number of parameters

# Print the results
print("Likelihood ratio test statistic:", lr_statistic)
print("p-value:", p_value)
if p_value < 0.05:
    print("The null hypothesis is rejected at the 95% significance level.")
else:
    print("The null hypothesis is not rejected at the 95% significance level.")
```

    Likelihood ratio test statistic: 30.39002044699052
    p-value: 2.517044579761756e-07
    The null hypothesis is rejected at the 95% significance level.