numcodecs_asinh/
lib.rs

1//! [![CI Status]][workflow] [![MSRV]][repo] [![Latest Version]][crates.io] [![Rust Doc Crate]][docs.rs] [![Rust Doc Main]][docs]
2//!
3//! [CI Status]: https://img.shields.io/github/actions/workflow/status/juntyr/numcodecs-rs/ci.yml?branch=main
4//! [workflow]: https://github.com/juntyr/numcodecs-rs/actions/workflows/ci.yml?query=branch%3Amain
5//!
6//! [MSRV]: https://img.shields.io/badge/MSRV-1.82.0-blue
7//! [repo]: https://github.com/juntyr/numcodecs-rs
8//!
9//! [Latest Version]: https://img.shields.io/crates/v/numcodecs-asinh
10//! [crates.io]: https://crates.io/crates/numcodecs-asinh
11//!
12//! [Rust Doc Crate]: https://img.shields.io/docsrs/numcodecs-asinh
13//! [docs.rs]: https://docs.rs/numcodecs-asinh/
14//!
15//! [Rust Doc Main]: https://img.shields.io/badge/docs-main-blue
16//! [docs]: https://juntyr.github.io/numcodecs-rs/numcodecs_asinh
17//!
18//! `$\text{asinh(x)}$` codec implementation for the [`numcodecs`] API.
19
20use ndarray::{Array, ArrayBase, ArrayView, ArrayViewMut, Data, Dimension, Zip};
21use num_traits::{Float, Signed};
22use numcodecs::{
23    AnyArray, AnyArrayAssignError, AnyArrayDType, AnyArrayView, AnyArrayViewMut, AnyCowArray,
24    Codec, StaticCodec, StaticCodecConfig, StaticCodecVersion,
25};
26use schemars::JsonSchema;
27use serde::{Deserialize, Serialize};
28use thiserror::Error;
29
30#[derive(Clone, Serialize, Deserialize, JsonSchema)]
31#[serde(deny_unknown_fields)]
32/// Asinh codec, which applies a quasi-logarithmic transformation on encoding.
33///
34/// For values close to zero that are within the codec's `linear_width`, the
35/// transform is close to linear. For values of larger magnitudes, the
36/// transform is asymptotically logarithmic. Unlike a logarithmic transform,
37/// this codec supports all finite values, including negative values and zero.
38///
39/// In detail, the codec calculates
40/// `$c = w \cdot \text{asinh}\left( \frac{x}{w} \right)$`
41/// on encoding and
42/// `$d = w \cdot \text{sinh}\left( \frac{c}{w} \right)$`
43/// on decoding, where `$w$` is the codec's `linear_width`.
44///
45/// The codec only supports finite floating point numbers.
46pub struct AsinhCodec {
47    /// The width of the close-to-zero input value range where the transform is
48    /// nearly linear
49    pub linear_width: f64,
50    /// The codec's encoding format version. Do not provide this parameter explicitly.
51    #[serde(default, rename = "_version")]
52    pub version: StaticCodecVersion<1, 0, 0>,
53}
54
55impl Codec for AsinhCodec {
56    type Error = AsinhCodecError;
57
58    fn encode(&self, data: AnyCowArray) -> Result<AnyArray, Self::Error> {
59        match data {
60            #[expect(clippy::cast_possible_truncation)]
61            AnyCowArray::F32(data) => Ok(AnyArray::F32(asinh(data, self.linear_width as f32)?)),
62            AnyCowArray::F64(data) => Ok(AnyArray::F64(asinh(data, self.linear_width)?)),
63            encoded => Err(AsinhCodecError::UnsupportedDtype(encoded.dtype())),
64        }
65    }
66
67    fn decode(&self, encoded: AnyCowArray) -> Result<AnyArray, Self::Error> {
68        match encoded {
69            #[expect(clippy::cast_possible_truncation)]
70            AnyCowArray::F32(encoded) => {
71                Ok(AnyArray::F32(sinh(encoded, self.linear_width as f32)?))
72            }
73            AnyCowArray::F64(encoded) => Ok(AnyArray::F64(sinh(encoded, self.linear_width)?)),
74            encoded => Err(AsinhCodecError::UnsupportedDtype(encoded.dtype())),
75        }
76    }
77
78    fn decode_into(
79        &self,
80        encoded: AnyArrayView,
81        decoded: AnyArrayViewMut,
82    ) -> Result<(), Self::Error> {
83        match (encoded, decoded) {
84            #[expect(clippy::cast_possible_truncation)]
85            (AnyArrayView::F32(encoded), AnyArrayViewMut::F32(decoded)) => {
86                sinh_into(encoded, decoded, self.linear_width as f32)
87            }
88            (AnyArrayView::F64(encoded), AnyArrayViewMut::F64(decoded)) => {
89                sinh_into(encoded, decoded, self.linear_width)
90            }
91            (encoded @ (AnyArrayView::F32(_) | AnyArrayView::F64(_)), decoded) => {
92                Err(AsinhCodecError::MismatchedDecodeIntoArray {
93                    source: AnyArrayAssignError::DTypeMismatch {
94                        src: encoded.dtype(),
95                        dst: decoded.dtype(),
96                    },
97                })
98            }
99            (encoded, _decoded) => Err(AsinhCodecError::UnsupportedDtype(encoded.dtype())),
100        }
101    }
102}
103
104impl StaticCodec for AsinhCodec {
105    const CODEC_ID: &'static str = "asinh.rs";
106
107    type Config<'de> = Self;
108
109    fn from_config(config: Self::Config<'_>) -> Self {
110        config
111    }
112
113    fn get_config(&self) -> StaticCodecConfig<Self> {
114        StaticCodecConfig::from(self)
115    }
116}
117
118#[derive(Debug, Error)]
119/// Errors that may occur when applying the [`AsinhCodec`].
120pub enum AsinhCodecError {
121    /// [`AsinhCodec`] does not support the dtype
122    #[error("Asinh does not support the dtype {0}")]
123    UnsupportedDtype(AnyArrayDType),
124    /// [`AsinhCodec`] does not support non-finite (infinite or NaN) floating
125    /// point data
126    #[error("Asinh does not support non-finite (infinite or NaN) floating point data")]
127    NonFiniteData,
128    /// [`AsinhCodec`] cannot decode into the provided array
129    #[error("Asinh cannot decode into the provided array")]
130    MismatchedDecodeIntoArray {
131        /// The source of the error
132        #[from]
133        source: AnyArrayAssignError,
134    },
135}
136
137/// Compute `$w \cdot \text{asinh}\left( \frac{x}{w} \right)$` over the
138/// elements of the input `data` array.
139///
140/// # Errors
141///
142/// Errors with
143/// - [`AsinhCodecError::NonFiniteData`] if any data element is non-finite
144///   (infinite or NaN)
145pub fn asinh<T: Float + Signed, S: Data<Elem = T>, D: Dimension>(
146    data: ArrayBase<S, D>,
147    linear_width: T,
148) -> Result<Array<T, D>, AsinhCodecError> {
149    if !Zip::from(&data).all(|x| x.is_finite()) {
150        return Err(AsinhCodecError::NonFiniteData);
151    }
152
153    let mut data = data.into_owned();
154    data.mapv_inplace(|x| (x / linear_width).asinh() * linear_width);
155
156    Ok(data)
157}
158
159/// Compute `$w \cdot \text{sinh}\left( \frac{x}{w} \right)$` over the
160/// elements of the input `data` array.
161///
162/// # Errors
163///
164/// Errors with
165/// - [`AsinhCodecError::NonFiniteData`] if any data element is non-finite
166///   (infinite or NaN)
167pub fn sinh<T: Float, S: Data<Elem = T>, D: Dimension>(
168    data: ArrayBase<S, D>,
169    linear_width: T,
170) -> Result<Array<T, D>, AsinhCodecError> {
171    if !Zip::from(&data).all(|x| x.is_finite()) {
172        return Err(AsinhCodecError::NonFiniteData);
173    }
174
175    let mut data = data.into_owned();
176    data.mapv_inplace(|x| (x / linear_width).sinh() * linear_width);
177
178    Ok(data)
179}
180
181#[expect(clippy::needless_pass_by_value)]
182/// Compute `$w \cdot \text{sinh}\left( \frac{x}{w} \right)$` over the
183/// elements of the input `data` array and write them into the `out`put array.
184///
185/// # Errors
186///
187/// Errors with
188/// - [`AsinhCodecError::NonFiniteData`] if any data element is non-finite
189///   (infinite or NaN)
190/// - [`AsinhCodecError::MismatchedDecodeIntoArray`] if the `data` array's shape
191///   does not match the `out`put array's shape
192pub fn sinh_into<T: Float, D: Dimension>(
193    data: ArrayView<T, D>,
194    mut out: ArrayViewMut<T, D>,
195    linear_width: T,
196) -> Result<(), AsinhCodecError> {
197    if data.shape() != out.shape() {
198        return Err(AsinhCodecError::MismatchedDecodeIntoArray {
199            source: AnyArrayAssignError::ShapeMismatch {
200                src: data.shape().to_vec(),
201                dst: out.shape().to_vec(),
202            },
203        });
204    }
205
206    if !Zip::from(&data).all(|x| x.is_finite()) {
207        return Err(AsinhCodecError::NonFiniteData);
208    }
209
210    // iteration must occur in synchronised (standard) order
211    for (d, o) in data.iter().zip(out.iter_mut()) {
212        *o = ((*d) / linear_width).sinh() * linear_width;
213    }
214
215    Ok(())
216}
217
218#[cfg(test)]
219mod tests {
220    use super::*;
221
222    #[test]
223    fn roundtrip() -> Result<(), AsinhCodecError> {
224        let data = (-1000..1000).map(f64::from).collect::<Vec<_>>();
225        let data = Array::from_vec(data);
226
227        let encoded = asinh(data.view(), 1.0)?;
228
229        for (r, e) in data.iter().zip(encoded.iter()) {
230            assert_eq!((*r).asinh().to_bits(), (*e).to_bits());
231        }
232
233        let decoded = sinh(encoded, 1.0)?;
234
235        for (r, d) in data.iter().zip(decoded.iter()) {
236            assert!(((*r) - (*d)).abs() < 1e-12);
237        }
238
239        Ok(())
240    }
241
242    #[test]
243    fn roundtrip_widths() -> Result<(), AsinhCodecError> {
244        let data = (-1000..1000).map(f64::from).collect::<Vec<_>>();
245        let data = Array::from_vec(data);
246
247        for linear_width in [-100.0, -10.0, -1.0, -0.1, 0.1, 1.0, 10.0, 100.0] {
248            let encoded = asinh(data.view(), linear_width)?;
249            let decoded = sinh(encoded, linear_width)?;
250
251            for (r, d) in data.iter().zip(decoded.iter()) {
252                assert!(((*r) - (*d)).abs() < 1e-12);
253            }
254        }
255
256        Ok(())
257    }
258}