From 431b410441fc06b508d862030ff27709335bd96c Mon Sep 17 00:00:00 2001 From: xuewc Date: Sat, 22 Feb 2025 00:40:34 +0800 Subject: [PATCH 1/4] fix: handle sparse matrix reading for big-endian byte order --- src/elisa/data/base.py | 10 ++++++++-- src/elisa/data/ogip.py | 6 ++++++ tests/test_data.py | 27 +++++++++++++++++++++++++++ 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/src/elisa/data/base.py b/src/elisa/data/base.py index 98fa0cf3..3da33991 100644 --- a/src/elisa/data/base.py +++ b/src/elisa/data/base.py @@ -1138,8 +1138,14 @@ def __init__( photon_egrid = np.array(photon_egrid, dtype=np.float64, order='C') channel_emin = np.array(channel_emin, dtype=np.float64, order='C') channel_emax = np.array(channel_emax, dtype=np.float64, order='C') - response_matrix = coo_array(response_matrix) - response_matrix.eliminate_zeros() + # if the byte order is not native, + # then convert it to be compatible with scipy.sparse + if not isinstance(response_matrix, sparray): + if response_matrix.dtype.byteorder != '=': + order = response_matrix.dtype.newbyteorder() + response_matrix = response_matrix.astype(order) + response_matrix = coo_array(response_matrix) + response_matrix.eliminate_zeros() channel = np.array(channel, dtype=str, order='C') photon_egrid_shape = np.shape(photon_egrid) diff --git a/src/elisa/data/ogip.py b/src/elisa/data/ogip.py index 3f380056..b580ec8d 100644 --- a/src/elisa/data/ogip.py +++ b/src/elisa/data/ogip.py @@ -629,6 +629,12 @@ def _read_response(self, file: str, response_id: int) -> dict: else: reduced_matrix = np.hstack(matrix) + # if the byte order is not native, + # then convert it to be compatible with scipy.sparse + if reduced_matrix.dtype.byteorder != '=': + order = reduced_matrix.dtype.newbyteorder() + reduced_matrix = reduced_matrix.astype(order) + sparse_matrix = coo_array( (reduced_matrix, (rows, cols)), shape=(len(response_data), channel_number), diff --git a/tests/test_data.py b/tests/test_data.py index 46159070..6ae3090b 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -1,6 +1,8 @@ import numpy as np +import pytest from elisa.data.grouping import significance_gv, significance_lima +from elisa.data.ogip import Response, ResponseData from elisa.models import PowerLaw @@ -117,3 +119,28 @@ def test_data_grouping(): data.spec_counts, data.back_counts, data.back_errors, data.back_ratio ) assert np.all(sig >= scale) + + +@pytest.mark.parametrize( + 'file', + [ + '../docs/notebooks/data/P011160500104_LE.rsp', + '../docs/notebooks/data/P011160500104_ME.rsp', + '../docs/notebooks/data/P011160500104_HE.rsp', + ], +) +def test_load_response(file): + # Test Response against big-endian files + assert np.all(Response(file).channel_fwhm > 0) + + +def test_response(): + # Test ResponseData against different endianness + photon_egrid = np.linspace(1.0, 100.0, 101) + channel = np.arange(100) + channel_emin = photon_egrid[:-1] + channel_emax = photon_egrid[1:] + matrix1 = np.eye(100).astype('f4') + ResponseData(photon_egrid, channel_emin, channel_emax, matrix1, channel) + ResponseData(photon_egrid, channel_emin, channel_emax, matrix2, channel) From 1141ad9ad51f5e0e96edf2ec6c3f0c6398de18c7 Mon Sep 17 00:00:00 2001 From: xuewc Date: Sat, 22 Feb 2025 00:55:40 +0800 Subject: [PATCH 2/4] Improve code --- src/elisa/data/base.py | 8 +++----- src/elisa/data/ogip.py | 8 +++----- src/elisa/util/misc.py | 10 ++++++++++ tests/test_data.py | 23 ++++++++++++++++------- 4 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/elisa/data/base.py b/src/elisa/data/base.py index 3da33991..cd12ed0f 100644 --- a/src/elisa/data/base.py +++ b/src/elisa/data/base.py @@ -21,6 +21,7 @@ significance_lima, ) from elisa.plot.misc import get_colors +from elisa.util.misc import to_native_byteorder if TYPE_CHECKING: NDArray = np.ndarray @@ -1138,12 +1139,9 @@ def __init__( photon_egrid = np.array(photon_egrid, dtype=np.float64, order='C') channel_emin = np.array(channel_emin, dtype=np.float64, order='C') channel_emax = np.array(channel_emax, dtype=np.float64, order='C') - # if the byte order is not native, - # then convert it to be compatible with scipy.sparse if not isinstance(response_matrix, sparray): - if response_matrix.dtype.byteorder != '=': - order = response_matrix.dtype.newbyteorder() - response_matrix = response_matrix.astype(order) + # convert to native byteorder to be compatible with scipy.sparse + response_matrix = to_native_byteorder(response_matrix) response_matrix = coo_array(response_matrix) response_matrix.eliminate_zeros() channel = np.array(channel, dtype=str, order='C') diff --git a/src/elisa/data/ogip.py b/src/elisa/data/ogip.py index b580ec8d..836367f6 100644 --- a/src/elisa/data/ogip.py +++ b/src/elisa/data/ogip.py @@ -9,6 +9,7 @@ from scipy.sparse import coo_array from elisa.data.base import ObservationData, ResponseData, SpectrumData +from elisa.util.misc import to_native_byteorder if TYPE_CHECKING: NDArray = np.ndarray @@ -629,11 +630,8 @@ def _read_response(self, file: str, response_id: int) -> dict: else: reduced_matrix = np.hstack(matrix) - # if the byte order is not native, - # then convert it to be compatible with scipy.sparse - if reduced_matrix.dtype.byteorder != '=': - order = reduced_matrix.dtype.newbyteorder() - reduced_matrix = reduced_matrix.astype(order) + # convert to native byteorder to be compatible with scipy.sparse + reduced_matrix = to_native_byteorder(reduced_matrix) sparse_matrix = coo_array( (reduced_matrix, (rows, cols)), diff --git a/src/elisa/util/misc.py b/src/elisa/util/misc.py index b78a6791..71a81227 100644 --- a/src/elisa/util/misc.py +++ b/src/elisa/util/misc.py @@ -22,6 +22,8 @@ from collections.abc import Sequence from typing import Callable, Literal, TypeVar + from numpy import ndarray as NDArray + from elisa.util.typing import CompEval T = TypeVar('T') @@ -571,3 +573,11 @@ def _wrapper_progress_bar(i, vals): return _wrapper_progress_bar return progress_bar_fori_loop + + +def to_native_byteorder(arr: NDArray) -> NDArray: + """Convert an array to native byte order.""" + if arr.dtype.byteorder != '=': + return arr.astype(arr.dtype.newbyteorder('=')) + else: + return arr diff --git a/tests/test_data.py b/tests/test_data.py index 6ae3090b..7c7ced69 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from jax.experimental.sparse import BCSR from elisa.data.grouping import significance_gv, significance_lima from elisa.data.ogip import Response, ResponseData @@ -130,17 +131,25 @@ def test_data_grouping(): ], ) def test_load_response(file): - # Test Response against big-endian files - assert np.all(Response(file).channel_fwhm > 0) + # test Response against big-endian files + rsp = Response(file) + # test if the response matrix can be converted to a BCSR matrix in JAX + BCSR.from_scipy_sparse(rsp.matrix) + assert np.all(rsp.channel_fwhm > 0) def test_response(): - # Test ResponseData against different endianness + # test ResponseData against different endianness photon_egrid = np.linspace(1.0, 100.0, 101) channel = np.arange(100) channel_emin = photon_egrid[:-1] channel_emax = photon_egrid[1:] - matrix1 = np.eye(100).astype('f4') - ResponseData(photon_egrid, channel_emin, channel_emax, matrix1, channel) - ResponseData(photon_egrid, channel_emin, channel_emax, matrix2, channel) + mat1 = np.eye(100).astype('f4') + r1 = ResponseData(photon_egrid, channel_emin, channel_emax, mat1, channel) + r2 = ResponseData(photon_egrid, channel_emin, channel_emax, mat2, channel) + # test if the response matrix can be converted to a BCSR matrix in JAX + BCSR.from_scipy_sparse(r1.matrix) + BCSR.from_scipy_sparse(r2.matrix) + assert np.all(r1.matrix == r2.matrix) + assert np.all(r1.channel_fwhm == r2.channel_fwhm) From a732653e0a9c1a936fe754d54c31252b3f7b806c Mon Sep 17 00:00:00 2001 From: xuewc Date: Sat, 22 Feb 2025 01:09:22 +0800 Subject: [PATCH 3/4] Fix test --- tests/test_data.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/tests/test_data.py b/tests/test_data.py index 7c7ced69..c0babab9 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -134,8 +134,11 @@ def test_load_response(file): # test Response against big-endian files rsp = Response(file) # test if the response matrix can be converted to a BCSR matrix in JAX - BCSR.from_scipy_sparse(rsp.matrix) assert np.all(rsp.channel_fwhm > 0) + assert np.array_equal( + rsp.matrix, + BCSR.from_scipy_sparse(rsp.sparse_matrix).todense(), + ) def test_response(): @@ -149,7 +152,9 @@ def test_response(): r1 = ResponseData(photon_egrid, channel_emin, channel_emax, mat1, channel) r2 = ResponseData(photon_egrid, channel_emin, channel_emax, mat2, channel) # test if the response matrix can be converted to a BCSR matrix in JAX - BCSR.from_scipy_sparse(r1.matrix) - BCSR.from_scipy_sparse(r2.matrix) - assert np.all(r1.matrix == r2.matrix) + for r in [r1, r2]: + assert np.array_equal( + r.matrix, BCSR.from_scipy_sparse(r.sparse_matrix).todense() + ) + assert np.all(r1.channel == r2.channel) assert np.all(r1.channel_fwhm == r2.channel_fwhm) From 17de58ded870ea5ce7a98df2b1eaca0b67c091cc Mon Sep 17 00:00:00 2001 From: xuewc Date: Sat, 22 Feb 2025 01:15:27 +0800 Subject: [PATCH 4/4] Fix test --- tests/test_data.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_data.py b/tests/test_data.py index c0babab9..15dfc822 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -125,9 +125,9 @@ def test_data_grouping(): @pytest.mark.parametrize( 'file', [ - '../docs/notebooks/data/P011160500104_LE.rsp', - '../docs/notebooks/data/P011160500104_ME.rsp', - '../docs/notebooks/data/P011160500104_HE.rsp', + 'docs/notebooks/data/P011160500104_LE.rsp', + 'docs/notebooks/data/P011160500104_ME.rsp', + 'docs/notebooks/data/P011160500104_HE.rsp', ], ) def test_load_response(file):