conspire/math/tensor/quantity/
mod.rs

1#[cfg(test)]
2mod test;
3
4pub(crate) mod sparse_vec;
5pub(crate) mod sparse_vec_2d;
6pub(crate) mod vec;
7
8use super::{
9    Differentiable, Erase, Hessian, HessianBlock, Jacobian, Solution, SquareMatrix, Tensor,
10    TensorArray, Vector, rank_0::TensorRank0,
11};
12use crate::math::{TensorList, assert::FiniteDifference};
13use crate::units::{Dimensionless, UnitDiv, UnitHalves, UnitInv, UnitMul, UnitRoot};
14use std::{
15    cmp::Ordering,
16    fmt::{self, Display, Formatter},
17    marker::PhantomData,
18    ops::{Add, AddAssign, Div, DivAssign, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign},
19};
20
21/// Implemented only where the two types are the same, so that a unit may be
22/// named where it is discarded without allowing a different one.
23pub trait Is<T> {}
24
25impl<T> Is<T> for T {}
26
27/// A scalar carrying a physical unit.
28#[repr(transparent)]
29pub struct Quantity<U = Dimensionless>(TensorRank0, PhantomData<U>);
30
31impl<U> Quantity<U> {
32    /// Associated function for const type conversion.
33    pub(crate) const fn new(value: TensorRank0) -> Self {
34        Self(value, PhantomData)
35    }
36    /// Returns the value with its unit discarded.
37    pub const fn value(&self) -> TensorRank0 {
38        self.0
39    }
40    /// Returns the value, stating the unit being discarded.
41    pub const fn value_as<V>(&self) -> TensorRank0
42    where
43        U: Is<V>,
44    {
45        self.0
46    }
47}
48
49impl<U> Quantity<U> {
50    /// Returns the absolute value, which leaves the unit alone.
51    pub fn abs(self) -> Self {
52        Self::new(self.0.abs())
53    }
54    /// Returns whether the value is not a number.
55    pub fn is_nan(&self) -> bool {
56        self.0.is_nan()
57    }
58    /// Returns whether two quantities differ by more than the epsilon
59    /// relatively, at least one of them being large enough for that ratio to
60    /// mean anything.
61    pub fn differs(self, quantity: Self, epsilon: TensorRank0) -> bool {
62        ((self.0 / quantity.0 - 1.0).abs() >= epsilon
63            && (self.0.abs() >= epsilon || quantity.0.abs() >= epsilon))
64            || self.is_nan()
65            || quantity.is_nan()
66    }
67    /// Returns whether two quantities [differ](Self::differs) absolutely as
68    /// well as relatively.
69    pub fn differs_severely(self, quantity: Self, epsilon: TensorRank0) -> bool {
70        ((self.0 / quantity.0 - 1.0).abs() >= epsilon
71            && (self.0 - quantity.0).abs() >= epsilon
72            && (self.0.abs() >= epsilon || quantity.0.abs() >= epsilon))
73            || self.is_nan()
74            || quantity.is_nan()
75    }
76    /// Returns how two quantities of the same unit order, totally.
77    pub fn total_cmp(&self, quantity: &Self) -> Ordering {
78        self.0.total_cmp(&quantity.0)
79    }
80    /// Returns the lesser of two quantities of the same unit.
81    pub fn min(self, quantity: Self) -> Self {
82        Self::new(self.0.min(quantity.0))
83    }
84    /// Returns the greater of two quantities of the same unit.
85    pub fn max(self, quantity: Self) -> Self {
86        Self::new(self.0.max(quantity.0))
87    }
88}
89
90impl<U> Quantity<U>
91where
92    U: UnitHalves,
93{
94    /// Returns the quantity twice over, as the units the halves of a tuple take
95    /// from it.
96    pub const fn halves(
97        self,
98    ) -> (
99        Quantity<<U as UnitHalves>::First>,
100        Quantity<<U as UnitHalves>::Second>,
101    ) {
102        (Quantity::new(self.0), Quantity::new(self.0))
103    }
104}
105
106impl<U> Quantity<U>
107where
108    U: UnitRoot,
109{
110    /// Returns the square root, which carries the unit this one is the square of.
111    pub fn sqrt(self) -> Quantity<<U as UnitRoot>::Output> {
112        Quantity::new(self.0.sqrt())
113    }
114}
115
116impl Quantity<Dimensionless> {
117    /// Returns the smallest integer greater than or equal to the value.
118    pub fn ceil(self) -> Self {
119        Self::new(self.0.ceil())
120    }
121    /// Returns the largest integer less than or equal to the value.
122    pub fn floor(self) -> Self {
123        Self::new(self.0.floor())
124    }
125    /// Raises to an integer power.
126    pub fn powi(self, n: i32) -> Self {
127        Self::new(self.0.powi(n))
128    }
129    /// Raises to a power.
130    pub fn powf(self, n: TensorRank0) -> Self {
131        Self::new(self.0.powf(n))
132    }
133    /// Returns the natural logarithm.
134    pub fn ln(self) -> Self {
135        Self::new(self.0.ln())
136    }
137    /// Returns the base-2 logarithm.
138    pub fn log2(self) -> Self {
139        Self::new(self.0.log2())
140    }
141    /// Returns the exponential.
142    pub fn exp(self) -> Self {
143        Self::new(self.0.exp())
144    }
145    /// Returns the sine.
146    pub fn sin(self) -> Self {
147        Self::new(self.0.sin())
148    }
149    /// Returns the cosine.
150    pub fn cos(self) -> Self {
151        Self::new(self.0.cos())
152    }
153}
154
155impl<U> Clone for Quantity<U> {
156    fn clone(&self) -> Self {
157        *self
158    }
159}
160
161impl<U> Copy for Quantity<U> {}
162
163impl<U> fmt::Debug for Quantity<U> {
164    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
165        fmt::Debug::fmt(&self.0, f)
166    }
167}
168
169impl<U> Display for Quantity<U> {
170    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
171        Display::fmt(&self.0, f)
172    }
173}
174
175impl<U> PartialEq for Quantity<U> {
176    fn eq(&self, other: &Self) -> bool {
177        self.0 == other.0
178    }
179}
180
181impl<U> PartialOrd for Quantity<U> {
182    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
183        self.0.partial_cmp(&other.0)
184    }
185}
186
187impl<U> Neg for Quantity<U> {
188    type Output = Self;
189    fn neg(self) -> Self::Output {
190        Self::new(-self.0)
191    }
192}
193
194impl<U> Default for Quantity<U> {
195    fn default() -> Self {
196        Self::new(0.0)
197    }
198}
199
200impl<U> Add for Quantity<U> {
201    type Output = Self;
202    fn add(self, quantity: Self) -> Self::Output {
203        Self::new(self.0 + quantity.0)
204    }
205}
206
207impl<U> Add<&Self> for Quantity<U> {
208    type Output = Self;
209    fn add(self, quantity: &Self) -> Self::Output {
210        Self::new(self.0 + quantity.0)
211    }
212}
213
214impl<U> AddAssign<&Self> for Quantity<U> {
215    fn add_assign(&mut self, quantity: &Self) {
216        self.0 += quantity.0
217    }
218}
219
220impl<U> Sub<&Self> for Quantity<U> {
221    type Output = Self;
222    fn sub(self, quantity: &Self) -> Self::Output {
223        Self::new(self.0 - quantity.0)
224    }
225}
226
227impl<U> SubAssign<&Self> for Quantity<U> {
228    fn sub_assign(&mut self, quantity: &Self) {
229        self.0 -= quantity.0
230    }
231}
232
233impl<U> Add for &Quantity<U> {
234    type Output = Quantity<U>;
235    fn add(self, quantity: Self) -> Self::Output {
236        Quantity::new(self.0 + quantity.0)
237    }
238}
239
240impl<U> Add<Quantity<U>> for &Quantity<U> {
241    type Output = Quantity<U>;
242    fn add(self, quantity: Quantity<U>) -> Self::Output {
243        Quantity::new(self.0 + quantity.0)
244    }
245}
246
247impl<U> Sub<Quantity<U>> for &Quantity<U> {
248    type Output = Quantity<U>;
249    fn sub(self, quantity: Quantity<U>) -> Self::Output {
250        Quantity::new(self.0 - quantity.0)
251    }
252}
253
254impl<U> Div<TensorRank0> for &Quantity<U> {
255    type Output = Quantity<U>;
256    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
257        Quantity::new(self.0 / tensor_rank_0)
258    }
259}
260
261impl<U> Div<&TensorRank0> for &Quantity<U> {
262    type Output = Quantity<U>;
263    fn div(self, tensor_rank_0: &TensorRank0) -> Self::Output {
264        Quantity::new(self.0 / tensor_rank_0)
265    }
266}
267
268impl<U> Neg for &Quantity<U> {
269    type Output = Quantity<U>;
270    fn neg(self) -> Self::Output {
271        Quantity::new(-self.0)
272    }
273}
274
275impl<U> Sub for &Quantity<U> {
276    type Output = Quantity<U>;
277    fn sub(self, quantity: Self) -> Self::Output {
278        Quantity::new(self.0 - quantity.0)
279    }
280}
281
282impl<U> Mul<TensorRank0> for &Quantity<U> {
283    type Output = Quantity<U>;
284    fn mul(self, tensor_rank_0: TensorRank0) -> Self::Output {
285        Quantity::new(self.0 * tensor_rank_0)
286    }
287}
288
289impl<U> Mul<&TensorRank0> for &Quantity<U> {
290    type Output = Quantity<U>;
291    fn mul(self, tensor_rank_0: &TensorRank0) -> Self::Output {
292        Quantity::new(self.0 * tensor_rank_0)
293    }
294}
295
296impl<U> MulAssign<&TensorRank0> for Quantity<U> {
297    fn mul_assign(&mut self, tensor_rank_0: &TensorRank0) {
298        self.0 *= tensor_rank_0
299    }
300}
301
302impl<U> DivAssign<&TensorRank0> for Quantity<U> {
303    fn div_assign(&mut self, tensor_rank_0: &TensorRank0) {
304        self.0 /= tensor_rank_0
305    }
306}
307
308impl<U> AddAssign for Quantity<U> {
309    fn add_assign(&mut self, quantity: Self) {
310        self.0 += quantity.0
311    }
312}
313
314impl<U> Sub for Quantity<U> {
315    type Output = Self;
316    fn sub(self, quantity: Self) -> Self::Output {
317        Self::new(self.0 - quantity.0)
318    }
319}
320
321impl<U> SubAssign for Quantity<U> {
322    fn sub_assign(&mut self, quantity: Self) {
323        self.0 -= quantity.0
324    }
325}
326
327impl<U> Mul<TensorRank0> for Quantity<U> {
328    type Output = Self;
329    fn mul(self, tensor_rank_0: TensorRank0) -> Self::Output {
330        Self::new(self.0 * tensor_rank_0)
331    }
332}
333
334impl<U> Mul<&TensorRank0> for Quantity<U> {
335    type Output = Self;
336    fn mul(self, tensor_rank_0: &TensorRank0) -> Self::Output {
337        Self::new(self.0 * tensor_rank_0)
338    }
339}
340
341impl<U> MulAssign<TensorRank0> for Quantity<U> {
342    fn mul_assign(&mut self, tensor_rank_0: TensorRank0) {
343        self.0 *= tensor_rank_0
344    }
345}
346
347impl<U> Div<TensorRank0> for Quantity<U> {
348    type Output = Self;
349    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
350        Self::new(self.0 / tensor_rank_0)
351    }
352}
353
354impl<U> DivAssign<TensorRank0> for Quantity<U> {
355    fn div_assign(&mut self, tensor_rank_0: TensorRank0) {
356        self.0 /= tensor_rank_0
357    }
358}
359
360impl<U> Mul<Quantity<U>> for TensorRank0 {
361    type Output = Quantity<U>;
362    fn mul(self, quantity: Quantity<U>) -> Self::Output {
363        Quantity::new(self * quantity.0)
364    }
365}
366
367/// Scaling a bare scalar by a dimensionless quantity leaves it bare.
368impl Mul<Quantity<Dimensionless>> for &TensorRank0 {
369    type Output = TensorRank0;
370    fn mul(self, quantity: Quantity<Dimensionless>) -> Self::Output {
371        self * quantity.0
372    }
373}
374
375impl<U, V> Mul<Quantity<V>> for Quantity<U>
376where
377    U: UnitMul<V>,
378{
379    type Output = Quantity<<U as UnitMul<V>>::Output>;
380    fn mul(self, quantity: Quantity<V>) -> Self::Output {
381        Quantity::new(self.0 * quantity.0)
382    }
383}
384
385impl<U, V> Mul<Quantity<V>> for &Quantity<U>
386where
387    U: UnitMul<V>,
388{
389    type Output = Quantity<<U as UnitMul<V>>::Output>;
390    fn mul(self, quantity: Quantity<V>) -> Self::Output {
391        Quantity::new(self.0 * quantity.0)
392    }
393}
394
395impl<U, V> Div<Quantity<V>> for &Quantity<U>
396where
397    U: UnitDiv<V>,
398{
399    type Output = Quantity<<U as UnitDiv<V>>::Output>;
400    fn div(self, quantity: Quantity<V>) -> Self::Output {
401        Quantity::new(self.0 / quantity.value())
402    }
403}
404
405impl<U, V> Div<Quantity<V>> for Quantity<U>
406where
407    U: UnitDiv<V>,
408{
409    type Output = Quantity<<U as UnitDiv<V>>::Output>;
410    fn div(self, quantity: Quantity<V>) -> Self::Output {
411        Quantity::new(self.0 / quantity.0)
412    }
413}
414
415impl<U, V> Mul<&Quantity<V>> for Quantity<U>
416where
417    U: UnitMul<V>,
418{
419    type Output = Quantity<<U as UnitMul<V>>::Output>;
420    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
421        self * *quantity
422    }
423}
424
425impl<U, V> Mul<&Quantity<V>> for &Quantity<U>
426where
427    U: UnitMul<V>,
428{
429    type Output = Quantity<<U as UnitMul<V>>::Output>;
430    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
431        *self * *quantity
432    }
433}
434
435impl<V> Mul<&Quantity<V>> for TensorRank0 {
436    type Output = Quantity<V>;
437    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
438        Quantity::new(self * quantity.0)
439    }
440}
441
442impl<V> Mul<&Quantity<V>> for &TensorRank0 {
443    type Output = Quantity<V>;
444    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
445        Quantity::new(self * quantity.0)
446    }
447}
448
449impl Add<TensorRank0> for Quantity<Dimensionless> {
450    type Output = Self;
451    fn add(self, tensor_rank_0: TensorRank0) -> Self::Output {
452        Self::new(self.0 + tensor_rank_0)
453    }
454}
455
456impl Add<Quantity<Dimensionless>> for TensorRank0 {
457    type Output = Quantity<Dimensionless>;
458    fn add(self, quantity: Quantity<Dimensionless>) -> Self::Output {
459        Quantity::new(self + quantity.0)
460    }
461}
462
463impl Sub<TensorRank0> for Quantity<Dimensionless> {
464    type Output = Self;
465    fn sub(self, tensor_rank_0: TensorRank0) -> Self::Output {
466        Self::new(self.0 - tensor_rank_0)
467    }
468}
469
470impl Sub<Quantity<Dimensionless>> for TensorRank0 {
471    type Output = Quantity<Dimensionless>;
472    fn sub(self, quantity: Quantity<Dimensionless>) -> Self::Output {
473        Quantity::new(self - quantity.0)
474    }
475}
476
477impl PartialEq<TensorRank0> for Quantity<Dimensionless> {
478    fn eq(&self, tensor_rank_0: &TensorRank0) -> bool {
479        &self.0 == tensor_rank_0
480    }
481}
482
483impl PartialEq<Quantity<Dimensionless>> for TensorRank0 {
484    fn eq(&self, quantity: &Quantity<Dimensionless>) -> bool {
485        self == &quantity.0
486    }
487}
488
489impl PartialOrd<TensorRank0> for Quantity<Dimensionless> {
490    fn partial_cmp(&self, tensor_rank_0: &TensorRank0) -> Option<std::cmp::Ordering> {
491        self.0.partial_cmp(tensor_rank_0)
492    }
493}
494
495impl PartialOrd<Quantity<Dimensionless>> for TensorRank0 {
496    fn partial_cmp(&self, quantity: &Quantity<Dimensionless>) -> Option<std::cmp::Ordering> {
497        self.partial_cmp(&quantity.0)
498    }
499}
500
501impl<U> Div<Quantity<U>> for TensorRank0
502where
503    U: UnitInv,
504{
505    type Output = Quantity<<U as UnitInv>::Output>;
506    fn div(self, quantity: Quantity<U>) -> Self::Output {
507        Quantity::new(self / quantity.0)
508    }
509}
510
511impl<U> std::iter::Sum for Quantity<U> {
512    fn sum<I>(iter: I) -> Self
513    where
514        I: Iterator<Item = Self>,
515    {
516        Self::new(iter.map(|quantity| quantity.0).sum())
517    }
518}
519
520impl<'a, U> std::iter::Sum<&'a Quantity<U>> for Quantity<U> {
521    fn sum<I>(iter: I) -> Self
522    where
523        I: Iterator<Item = &'a Quantity<U>>,
524    {
525        Self::new(iter.map(|quantity| quantity.0).sum())
526    }
527}
528
529impl<U> Erase for Quantity<U> {
530    type Erased = TensorRank0;
531    fn erase(&self) -> &Self::Erased {
532        &self.0
533    }
534}
535
536impl<U> Tensor for Quantity<U> {
537    type Item = Self;
538    type Unit = U;
539    fn error_count_zero(&self, tol_abs: TensorRank0, tol_rel: TensorRank0) -> Option<usize> {
540        self.0.error_count_zero(tol_abs, tol_rel)
541    }
542    fn error_count(
543        &self,
544        other: &Self,
545        tol_abs: TensorRank0,
546        tol_rel: TensorRank0,
547    ) -> Option<usize> {
548        self.0.error_count(&other.0, tol_abs, tol_rel)
549    }
550    fn full_contraction(&self, quantity: &Self) -> TensorRank0 {
551        self.0 * quantity.0
552    }
553    fn is_zero(&self) -> bool {
554        self.0 == 0.0
555    }
556    fn iter(&self) -> impl Iterator<Item = &Self::Item> {
557        std::slice::from_ref(self).iter()
558    }
559    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item> {
560        std::slice::from_mut(self).iter_mut()
561    }
562    fn len(&self) -> usize {
563        1
564    }
565    fn norm_inf(&self) -> Quantity<U> {
566        Self::new(self.0.abs())
567    }
568    fn norm_l1(&self) -> Quantity<U> {
569        Self::new(self.0.abs())
570    }
571    fn norm_p_sum(&self, p: TensorRank0) -> TensorRank0 {
572        self.0.abs().powf(p)
573    }
574    fn size(&self) -> usize {
575        1
576    }
577    fn sub_abs(&self, other: &Self) -> Self {
578        Self::new((self.0 - other.0).abs())
579    }
580    fn sub_rel(&self, other: &Self) -> Self {
581        Self::new(self.0.sub_rel(&other.0))
582    }
583}
584
585impl<U> TensorArray for Quantity<U> {
586    type Array = TensorRank0;
587    type Item = Self;
588    fn as_array(&self) -> Self::Array {
589        self.0
590    }
591    fn identity() -> Self {
592        Self::new(1.0)
593    }
594    fn zero() -> Self {
595        Self::new(0.0)
596    }
597}
598
599impl<U> FiniteDifference for Quantity<U> {
600    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
601        self.0.error_fd(&comparator.0, epsilon)
602    }
603}
604
605impl<U, const N: usize> FiniteDifference for TensorList<Quantity<U>, N> {
606    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
607        error_fd_over(self.iter().zip(comparator.iter()), epsilon)
608    }
609}
610
611impl<U, const M: usize, const N: usize> FiniteDifference
612    for TensorList<TensorList<Quantity<U>, N>, M>
613{
614    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
615        error_fd_over(
616            self.iter()
617                .zip(comparator.iter())
618                .flat_map(|(entry, comparator_entry)| entry.iter().zip(comparator_entry.iter())),
619            epsilon,
620        )
621    }
622}
623
624fn error_fd_over<'a, U: 'a>(
625    entries: impl Iterator<Item = (&'a Quantity<U>, &'a Quantity<U>)>,
626    epsilon: TensorRank0,
627) -> Option<(bool, usize)> {
628    let error_count = entries
629        .filter_map(|(entry, comparator_entry)| entry.error_fd(comparator_entry, epsilon))
630        .map(|(_, count)| count)
631        .sum();
632    if error_count > 0 {
633        Some((true, error_count))
634    } else {
635        None
636    }
637}
638
639impl<U> Solution for Quantity<U> {
640    fn decrement_from(&mut self, other: &Vector) {
641        self.0 -= other[0]
642    }
643    fn decrement_from_chained(&mut self, other: &mut Vector, vector: &Vector) {
644        self.0 -= vector[0];
645        other
646            .iter_mut()
647            .zip(vector.iter().skip(1))
648            .for_each(|(entry_i, vector_i)| *entry_i -= vector_i)
649    }
650}
651
652impl<U> Hessian for Quantity<U> {
653    fn quadratic_form(&self, vector: &Vector) -> TensorRank0 {
654        self.0 * vector[0] * vector[0]
655    }
656    fn entry(&self, _row: usize, _column: usize) -> TensorRank0 {
657        unimplemented!()
658    }
659    fn fill_into(self, _square_matrix: &mut SquareMatrix) {
660        unimplemented!()
661    }
662}
663
664/// A quantity is a 1x1 block.
665impl<U> HessianBlock for Quantity<U> {
666    fn entry(&self, _row: usize, _column: usize) -> TensorRank0 {
667        self.0
668    }
669    fn height(&self) -> usize {
670        1
671    }
672    fn width(&self) -> usize {
673        1
674    }
675    fn fill_into_block<M>(&self, matrix: &mut M, row: usize, column: usize)
676    where
677        M: IndexMut<usize, Output = Vector>,
678    {
679        matrix[row][column] = self.0
680    }
681}
682
683impl<U> Jacobian for Quantity<U> {
684    fn fill_into(&self, vector: &mut Vector) {
685        vector[0] = self.0
686    }
687    fn fill_into_chained(self, other: Vector, vector: &mut Vector) {
688        vector[0] = self.0;
689        other
690            .into_iter()
691            .zip(vector.iter_mut().skip(1))
692            .for_each(|(entry_i, vector_i)| *vector_i = entry_i)
693    }
694}
695
696impl<U> Sub<Vector> for Quantity<U> {
697    type Output = Self;
698    fn sub(self, vector: Vector) -> Self::Output {
699        Self::new(self.0 - vector[0])
700    }
701}
702
703impl<U> Sub<&Vector> for Quantity<U> {
704    type Output = Self;
705    fn sub(self, vector: &Vector) -> Self::Output {
706        Self::new(self.0 - vector[0])
707    }
708}
709
710impl<U> From<Vector> for Quantity<U> {
711    fn from(vector: Vector) -> Self {
712        Self::new(vector[0])
713    }
714}
715
716impl<U, T> Differentiable<T> for Quantity<U>
717where
718    U: UnitDiv<T>,
719{
720    type Derivative = Quantity<<U as UnitDiv<T>>::Output>;
721}
conspire/math/tensor/quantity/mod.rs

conspire/math/tensor/quantity/
mod.rs
