RadClass/tests/test_DiffSpectra.py at 3b022a0eaf756e0cd7db1e994def4bcd0a536b2f · stompsjo/RadClass · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
import numpy as np
import pytest
import os
from datetime import datetime, timedelta

from RadClass.Processor import Processor
from RadClass.DiffSpectra import DiffSpectra
import tests.test_data as test_data

# initialize sample data
start_date = datetime(2019, 2, 2)
delta = timedelta(seconds=1)
timestamps = np.arange(start_date,
                       start_date + (test_data.timesteps * delta),
                       delta).astype('datetime64[s]').astype('float64')

live = np.full((len(timestamps),), test_data.livetime)
spectra = np.arange(test_data.timesteps)
spectra = np.full((test_data.energy_bins, spectra.shape[0]), spectra).T


@pytest.fixture(scope="module", autouse=True)
def init_test_file():
    # create sample test file with above simulated data
    yield test_data.create_file(live, timestamps, spectra)
    os.remove(test_data.filename)


def test_difference():
    stride = 10
    integration = 10

    # run handler script with analysis parameter
    # small stride since there are less data samples in test_data
    diff_stride = 2
    post_analysis = DiffSpectra(stride=diff_stride)
    classifier = Processor(stride, integration, test_data.datapath,
                           test_data.filename, post_analysis=post_analysis,
                           store_data=True)
    classifier.run_all()

    diff_spectra = post_analysis.diff_spectra
    # DiffSpectra length test:
    # there should be one difference spectrum for each timestamp with
    # diff_stride number of spectra before it
    exp_len = classifier.storage.shape[0] - diff_stride
    np.testing.assert_equal(diff_spectra.shape[0], exp_len)

    # DiffSpectra value test:
    # for test_data, the minimum background will always be the first spectrum
    # in a window (because the spectra increase 1, 2, 3, 4, etc.)
    # therefore the diff spectra element values will always be
    # spectra[i] - spectra[i-diff_stride]
    # the counts in a given window are defined analytically for
    # the spectral structure in test_data as:
    # n1 = integration + diff_stride * stride
    # n2 = integration + (diff_stride - 1) * stride
    # n3 = integration
    # diff_value = (n1^2 + n1) - (n2^2 + n2) - (n3^2 + n3)
    # (the algebra is simplified below)
    diff_value = (2*diff_stride*stride**2) - stride**2 - integration**2 \
        + (2*stride*integration) + stride - integration
    exp_spectra = np.full((exp_len, test_data.energy_bins),
                          diff_value/(2*test_data.livetime))
    np.testing.assert_almost_equal(diff_spectra[:, 1:], exp_spectra, decimal=2)

    # DiffSpectra timestamp test:
    # there should be one difference spectrum for each timestamp with
    # diff_stride number of spectra before it (i.e. spectra for every timestamp
    # after the first diff_stride spectra)
    exp_ts = classifier.storage[diff_stride:, 0]
    np.testing.assert_equal(diff_spectra[:, 0], exp_ts)


def test_write():
    stride = 10
    integration = 10
    filename = 'DiffSpectra_test.csv'

    # run handler script with analysis parameter
    # small stride since there are less data samples in test_data
    diff_stride = 2
    post_analysis = DiffSpectra(stride=diff_stride)
    classifier = Processor(stride, integration, test_data.datapath,
                           test_data.filename, post_analysis=post_analysis,
                           store_data=True)
    classifier.run_all()
    post_analysis.write(filename)

    results = np.loadtxt(filename, delimiter=',')
    # 1 extra columns are required for timestamp
    # expected shape is only 1D because only 1 entry is expected
    obs = results.shape
    exp = (classifier.storage.shape[0] - diff_stride, test_data.energy_bins+1)
    np.testing.assert_equal(obs, exp)

    os.remove(filename)