Skip to content

Commit

Permalink
Merge pull request #142 from vuillaut/astropy_units
Browse files Browse the repository at this point in the history 
Use astropy units whenever possible
  • Loading branch information
@vuillaut
vuillaut authored Apr 12, 2021
2 parents feaed54 + acbe28e commit dd59f45
Show file tree
Hide file tree
Showing 9 changed files with 1,358 additions and 1,182 deletions.
704 changes: 430 additions & 274 deletions ctaplot/ana/ana.py
View file

Large diffs are not rendered by default.

167 changes: 125 additions & 42 deletions ctaplot/ana/tests/test_ana.py
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,21 @@
import ctaplot.ana.ana as ana
import numpy as np
import astropy.units as u

np.random.seed(42)


def test_logspace_decades_nbin():
ca = ana.logspace_decades_nbin(0.1, 10, n=9)
assert len(ca) == 19
assert ca[0] == 0.1
assert ca[-1] == 10
ca = ana.logspace_decades_nbin(0.1, 10., n=10)
assert len(ca) == 21
np.testing.assert_almost_equal(np.log10(ca), np.linspace(-1, 1, 21))
np.testing.assert_almost_equal(ca[0], 0.1)
np.testing.assert_almost_equal(ca[-1], 10)


def test_logbin_mean():
np.testing.assert_allclose(ana.logbin_mean(np.array([0.1, 1, 10])), [np.sqrt(0.1), np.sqrt(10)])
np.testing.assert_allclose(ana.logbin_mean(np.array([0.1, 1, 10])*u.m), np.array([np.sqrt(0.1), np.sqrt(10)])*u.m)


def test_class_cta_requirement():
Expand All @@ -30,13 +37,13 @@ def test_class_cta_performance():


def test_impact_resolution_per_energy():
true_x = np.random.rand(100) * 1000
true_y = np.random.rand(100) * 1000
reco_x = true_x + 1
reco_y = true_y + 1
energy = np.logspace(-2, 2, 100)
E, R = ana.impact_resolution_per_energy(reco_x, reco_y, true_x, true_y, energy)
assert (np.isclose(R, np.sqrt(2))).all()
true_x = np.random.rand(100) * 1000 * u.m
true_y = np.random.rand(100) * 1000 * u.m
reco_x = true_x + 1 * u.m
reco_y = true_y + 1 * u.m
energy = np.logspace(-2, 2, 100) * u.TeV
E, R = ana.impact_resolution_per_energy(true_x, reco_x, true_y, reco_y, energy)
assert (np.isclose(R, np.sqrt(2) * u.m)).all()


def test_resolution():
Expand Down Expand Up @@ -99,7 +106,7 @@ def test_resolution_per_bin():

def test_resolution_per_bin_empty():
'''
testing for empty bins
testing for empty bins_x
'''
# For a normal distribution, the resolution at `percentile=68.27` is equal to 1 sigma
size = 1000000
Expand All @@ -122,20 +129,18 @@ def test_resolution_per_bin_empty():
def test_resolution_per_energy():
x = np.random.normal(size=100000, scale=1, loc=10)
y = 10 * np.ones(x.shape[0])
E = 10 ** (-3 + 6 * np.random.rand(x.shape[0]))
E = 10 ** (-3 + 6 * np.random.rand(x.shape[0])) * u.TeV
e_bin, res_e = ana.resolution_per_energy(y, x, E)
assert np.isclose(res_e, 1. / ana.relative_scaling(y, x).mean(), rtol=1e-1).all()


def test_power_law_integrated_distribution():
from ctaplot.ana.ana import power_law_integrated_distribution
emin = 50. # u.GeV
emax = 500.e3 # u.GeV
emin = 50. * u.GeV
emax = 500.e3 * u.GeV
Nevents = 1e6
spectral_index = -2.1
bins = np.logspace(np.log10(emin),
np.log10(emax),
40)
bins = np.geomspace(emin, emax, 40)

y = power_law_integrated_distribution(emin, emax, Nevents, spectral_index, bins)

Expand All @@ -150,15 +155,15 @@ def test_distance2d_resolution():
t = np.random.rand(size) * np.pi * 2
reco_x = 1 + 3 * np.cos(t)
reco_y = 1 + 3 * np.sin(t)
res, err_min, err_max = ana.distance2d_resolution(reco_x, reco_y, true_x, true_y,
res, err_min, err_max = ana.distance2d_resolution(true_x, reco_x, true_y, reco_y,
percentile=68.27, confidence_level=0.95, bias_correction=False)

np.testing.assert_equal(res, 3)

# with different bias on X and Y:
reco_x = 5 + 2 * np.cos(t)
reco_y = 7 + 2 * np.sin(t)
res, err_min, err_max = ana.distance2d_resolution(reco_x, reco_y, true_x, true_y,
res, err_min, err_max = ana.distance2d_resolution(true_x, reco_x, true_y, reco_y,
percentile=68.27, confidence_level=0.95, bias_correction=True)

assert np.isclose(res, 2, rtol=1e-1)
Expand All @@ -173,32 +178,32 @@ def test_distance2d_resolution_per_bin():
reco_x = 3 * np.cos(t)
reco_y = 3 * np.sin(t)

bin, res = ana.distance2d_resolution_per_bin(x, reco_x, reco_y, true_x, true_y, bins=10, bias_correction=True)
bin, res = ana.distance2d_resolution_per_bin(x, true_x, reco_x, true_y, reco_y, bins=10, bias_correction=True)

np.testing.assert_allclose(res, 3, rtol=1e-1)


def test_angular_resolution():
size = 10000
true_alt = np.random.rand(size)
true_az = np.random.rand(size)
true_alt = np.random.rand(size) * u.rad
true_az = np.random.rand(size) * u.rad
scale = 0.01
bias = 3

# test alt
reco_alt = true_alt + np.random.normal(loc=bias, scale=scale, size=size)
reco_alt = true_alt + np.random.normal(loc=bias, scale=scale, size=size)*u.rad
reco_az = true_az

assert np.isclose(ana.angular_resolution(reco_alt, reco_az, true_alt, true_az, bias_correction=True)[0],
scale,
assert np.isclose(ana.angular_resolution(true_alt, reco_alt, true_az, reco_az, bias_correction=True)[0],
scale * u.rad,
rtol=1e-1)

# test az
true_alt = np.zeros(size) # angular separation evolves linearly with az if alt=0
true_alt = np.zeros(size) * u.rad # angular separation evolves linearly with az if alt=0
reco_alt = true_alt
reco_az = true_az + np.random.normal(loc=-1, scale=scale, size=size)
assert np.isclose(ana.angular_resolution(reco_alt, reco_az, true_alt, true_az, bias_correction=True)[0],
scale,
reco_az = true_az + np.random.normal(loc=-1, scale=scale, size=size) * u.rad
assert np.isclose(ana.angular_resolution(true_alt, reco_alt, true_az, reco_az, bias_correction=True)[0],
scale * u.rad,
rtol=1e-1)


Expand All @@ -207,21 +212,23 @@ def test_angular_resolution_small_angles():
At small angles, the angular resolution should be equal to the distance2d resolution
"""
size = 1000
true_az = np.ones(size)
true_alt = np.random.rand(size)
reco_alt = true_alt + np.random.normal(1, 0.05, size)
true_az = np.ones(size) * u.rad
true_alt = np.random.rand(size) * u.rad
reco_alt = true_alt + np.random.normal(1, 0.05, size) * u.rad
reco_az = 2 * true_az
np.testing.assert_allclose(ana.angular_resolution(reco_alt, reco_az, true_alt, true_az, bias_correction=True),
ana.distance2d_resolution(reco_alt, reco_az, true_alt, true_az, bias_correction=True),
np.testing.assert_allclose(ana.angular_resolution(true_alt, reco_alt, true_az, reco_az, bias_correction=True),
ana.distance2d_resolution(true_alt, reco_alt, true_az, reco_az,
bias_correction=True),
rtol=1e-1,
)

true_az = np.random.rand(size)
true_alt = np.zeros(size)
true_az = np.random.rand(size) * u.rad
true_alt = np.zeros(size) * u.rad
reco_alt = true_alt
reco_az = true_az + np.random.normal(-3, 2, size)
np.testing.assert_allclose(ana.angular_resolution(reco_alt, reco_az, true_alt, true_az, bias_correction=True),
ana.distance2d_resolution(reco_alt, reco_az, true_alt, true_az, bias_correction=True),
reco_az = true_az + np.random.normal(-3, 2, size) * u.rad
np.testing.assert_allclose(ana.angular_resolution(true_alt, reco_alt, true_az, reco_az, bias_correction=True),
ana.distance2d_resolution(true_alt, reco_alt, true_az, reco_az,
bias_correction=True),
rtol=1e-1,
)

Expand All @@ -231,6 +238,12 @@ def test_bias_empty():
assert ana.bias(x, x) == 0


def test_bias():
reco = np.random.normal(loc=0, scale=0.00001, size=100) * u.TeV
true = np.random.normal(loc=1, scale=0.00001, size=100) * u.TeV
np.testing.assert_almost_equal(ana.bias(true, reco).to_value(u.TeV), -1, decimal=3)


def test_bias_per_bin():
size = 100000
true = np.ones(size)
Expand All @@ -244,7 +257,7 @@ def test_bias_per_energy():
size = 100000
true = np.ones(size)
reco = np.random.normal(loc=2, scale=0.5, size=size)
energy = np.logspace(-2, 2, size)
energy = np.logspace(-2, 2, size) * u.TeV
bins, bias = ana.bias_per_energy(true, reco, energy)
np.testing.assert_allclose(bias, 1, rtol=1e-1)

Expand All @@ -253,5 +266,75 @@ def test_get_magic_sensitivity():
from astropy.table.table import QTable
import astropy.units as u
table = ana.get_magic_sensitivity()
assert type(table) is QTable
assert isinstance(table, QTable)
assert table['e_min'][0] == 63 * u.GeV


def test_stat_per_energy():
size = 10000
energy = 10**(np.random.uniform(-2, 2, size=size)) * u.TeV
y = np.random.normal(loc=1, scale=0.01, size=size)
stat, edges, number = ana.stat_per_energy(energy, y, statistic='mean')
np.testing.assert_almost_equal(stat, 1, decimal=2)


def test_effective_area():
simu = 10**np.random.rand(100000) * u.TeV
reco = 10**np.random.rand(1000) * u.TeV
simu_area = 1*u.km**2
np.testing.assert_almost_equal(len(reco)/len(simu)*simu_area.to_value(u.m**2),
ana.effective_area(simu, reco, simu_area).to_value(u.m**2),
)


def test_effective_area_per_energy():
simu = np.random.uniform(1, 2, 1000000) * u.TeV
reco = np.random.uniform(1, 2, 10000) * u.TeV
simu_area = 1 * u.km ** 2
expected_area = len(reco) / len(simu) * simu_area.to_value(u.m ** 2)
bins, eff_area = ana.effective_area_per_energy(simu, reco, simu_area)
np.testing.assert_almost_equal(eff_area[eff_area!=0].to_value(u.m ** 2)/expected_area,
1,
decimal=2
)


def test_theta2():
true_alt = np.random.rand(10) * u.rad
reco_alt = np.random.rand(10) * u.rad
true_az = np.random.rand(10) * u.rad
reco_az = np.random.rand(10) * u.rad
t2 = ana.theta2(true_alt, reco_alt, true_az, reco_az)
assert t2.unit.is_equivalent(u.deg**2)


def test_energy_resolution_per_energy():
true_energy = 10**np.random.rand(50) * u.TeV
reco_energy = 10 ** np.random.rand(50) * u.TeV
ana.energy_resolution_per_energy(true_energy, reco_energy, bins=np.logspace(-1, 2)*u.TeV)


def test_energy_resolution():
true_energy = 10 ** np.random.rand(50) * u.TeV
reco_energy = 10 ** np.random.rand(50) * u.TeV
ana.energy_resolution(true_energy, reco_energy)


def test_angular_resolution_per_bin():
size=100
true_alt = np.random.rand(size) * u.rad
reco_alt = np.random.rand(size) * u.rad
true_az = np.random.rand(size) * u.rad
reco_az = np.random.rand(size) * u.rad
x = np.random.rand(size) * u.m
ana.angular_resolution_per_bin(true_alt, reco_alt, true_az, reco_az, x, bins=np.linspace(0, 1)*u.m)


def test_angular_resolution_per_energy():
size = 100
true_alt = np.random.rand(size) * u.rad
reco_alt = np.random.rand(size) * u.rad
true_az = np.random.rand(size) * u.rad
reco_az = np.random.rand(size) * u.rad
x = 10**np.random.rand(size) * u.TeV
ana.angular_resolution_per_energy(true_alt, reco_alt, true_az, reco_az, x, bins=np.logspace(0, 1) * u.TeV)
Loading

0 comments on commit dd59f45

Please sign in to comment.