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    Differentiate, Erase, Hessian, Jacobian, Solution, SquareMatrix, Tensor, TensorArray, Vector,
10    rank_0::TensorRank0,
11};
12use crate::math::{TensorList, assert::FiniteDifference};
13use crate::units::{Dimensionless, UnitDiv, UnitHalves, UnitInv, UnitMul};
14use std::{
15    cmp::Ordering,
16    fmt::{self, Display, Formatter},
17    marker::PhantomData,
18    ops::{Add, AddAssign, Div, DivAssign, 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 Quantity<Dimensionless> {
107    /// Returns the smallest integer greater than or equal to the value.
108    pub fn ceil(self) -> Self {
109        Self::new(self.0.ceil())
110    }
111    /// Returns the largest integer less than or equal to the value.
112    pub fn floor(self) -> Self {
113        Self::new(self.0.floor())
114    }
115    /// Raises to an integer power.
116    pub fn powi(self, n: i32) -> Self {
117        Self::new(self.0.powi(n))
118    }
119    /// Raises to a power.
120    pub fn powf(self, n: TensorRank0) -> Self {
121        Self::new(self.0.powf(n))
122    }
123    /// Returns the square root.
124    pub fn sqrt(self) -> Self {
125        Self::new(self.0.sqrt())
126    }
127    /// Returns the natural logarithm.
128    pub fn ln(self) -> Self {
129        Self::new(self.0.ln())
130    }
131    /// Returns the base-2 logarithm.
132    pub fn log2(self) -> Self {
133        Self::new(self.0.log2())
134    }
135    /// Returns the exponential.
136    pub fn exp(self) -> Self {
137        Self::new(self.0.exp())
138    }
139    /// Returns the sine.
140    pub fn sin(self) -> Self {
141        Self::new(self.0.sin())
142    }
143    /// Returns the cosine.
144    pub fn cos(self) -> Self {
145        Self::new(self.0.cos())
146    }
147}
148
149impl<U> Clone for Quantity<U> {
150    fn clone(&self) -> Self {
151        *self
152    }
153}
154
155impl<U> Copy for Quantity<U> {}
156
157impl<U> fmt::Debug for Quantity<U> {
158    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
159        fmt::Debug::fmt(&self.0, f)
160    }
161}
162
163impl<U> Display for Quantity<U> {
164    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
165        Display::fmt(&self.0, f)
166    }
167}
168
169impl<U> PartialEq for Quantity<U> {
170    fn eq(&self, other: &Self) -> bool {
171        self.0 == other.0
172    }
173}
174
175impl<U> PartialOrd for Quantity<U> {
176    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
177        self.0.partial_cmp(&other.0)
178    }
179}
180
181impl<U> Neg for Quantity<U> {
182    type Output = Self;
183    fn neg(self) -> Self::Output {
184        Self::new(-self.0)
185    }
186}
187
188impl<U> Default for Quantity<U> {
189    fn default() -> Self {
190        Self::new(0.0)
191    }
192}
193
194impl<U> Add for Quantity<U> {
195    type Output = Self;
196    fn add(self, quantity: Self) -> Self::Output {
197        Self::new(self.0 + quantity.0)
198    }
199}
200
201impl<U> Add<&Self> for Quantity<U> {
202    type Output = Self;
203    fn add(self, quantity: &Self) -> Self::Output {
204        Self::new(self.0 + quantity.0)
205    }
206}
207
208impl<U> AddAssign<&Self> for Quantity<U> {
209    fn add_assign(&mut self, quantity: &Self) {
210        self.0 += quantity.0
211    }
212}
213
214impl<U> Sub<&Self> for Quantity<U> {
215    type Output = Self;
216    fn sub(self, quantity: &Self) -> Self::Output {
217        Self::new(self.0 - quantity.0)
218    }
219}
220
221impl<U> SubAssign<&Self> for Quantity<U> {
222    fn sub_assign(&mut self, quantity: &Self) {
223        self.0 -= quantity.0
224    }
225}
226
227impl<U> Add for &Quantity<U> {
228    type Output = Quantity<U>;
229    fn add(self, quantity: Self) -> Self::Output {
230        Quantity::new(self.0 + quantity.0)
231    }
232}
233
234impl<U> Add<Quantity<U>> for &Quantity<U> {
235    type Output = Quantity<U>;
236    fn add(self, quantity: Quantity<U>) -> Self::Output {
237        Quantity::new(self.0 + quantity.0)
238    }
239}
240
241impl<U> Sub<Quantity<U>> for &Quantity<U> {
242    type Output = Quantity<U>;
243    fn sub(self, quantity: Quantity<U>) -> Self::Output {
244        Quantity::new(self.0 - quantity.0)
245    }
246}
247
248impl<U> Div<TensorRank0> for &Quantity<U> {
249    type Output = Quantity<U>;
250    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
251        Quantity::new(self.0 / tensor_rank_0)
252    }
253}
254
255impl<U> Div<&TensorRank0> for &Quantity<U> {
256    type Output = Quantity<U>;
257    fn div(self, tensor_rank_0: &TensorRank0) -> Self::Output {
258        Quantity::new(self.0 / tensor_rank_0)
259    }
260}
261
262impl<U> Neg for &Quantity<U> {
263    type Output = Quantity<U>;
264    fn neg(self) -> Self::Output {
265        Quantity::new(-self.0)
266    }
267}
268
269impl<U> Sub for &Quantity<U> {
270    type Output = Quantity<U>;
271    fn sub(self, quantity: Self) -> Self::Output {
272        Quantity::new(self.0 - quantity.0)
273    }
274}
275
276impl<U> Mul<TensorRank0> for &Quantity<U> {
277    type Output = Quantity<U>;
278    fn mul(self, tensor_rank_0: TensorRank0) -> Self::Output {
279        Quantity::new(self.0 * tensor_rank_0)
280    }
281}
282
283impl<U> Mul<&TensorRank0> for &Quantity<U> {
284    type Output = Quantity<U>;
285    fn mul(self, tensor_rank_0: &TensorRank0) -> Self::Output {
286        Quantity::new(self.0 * tensor_rank_0)
287    }
288}
289
290impl<U> MulAssign<&TensorRank0> for Quantity<U> {
291    fn mul_assign(&mut self, tensor_rank_0: &TensorRank0) {
292        self.0 *= tensor_rank_0
293    }
294}
295
296impl<U> DivAssign<&TensorRank0> for Quantity<U> {
297    fn div_assign(&mut self, tensor_rank_0: &TensorRank0) {
298        self.0 /= tensor_rank_0
299    }
300}
301
302impl<U> AddAssign for Quantity<U> {
303    fn add_assign(&mut self, quantity: Self) {
304        self.0 += quantity.0
305    }
306}
307
308impl<U> Sub for Quantity<U> {
309    type Output = Self;
310    fn sub(self, quantity: Self) -> Self::Output {
311        Self::new(self.0 - quantity.0)
312    }
313}
314
315impl<U> SubAssign for Quantity<U> {
316    fn sub_assign(&mut self, quantity: Self) {
317        self.0 -= quantity.0
318    }
319}
320
321impl<U> Mul<TensorRank0> for Quantity<U> {
322    type Output = Self;
323    fn mul(self, tensor_rank_0: TensorRank0) -> Self::Output {
324        Self::new(self.0 * tensor_rank_0)
325    }
326}
327
328impl<U> Mul<&TensorRank0> for Quantity<U> {
329    type Output = Self;
330    fn mul(self, tensor_rank_0: &TensorRank0) -> Self::Output {
331        Self::new(self.0 * tensor_rank_0)
332    }
333}
334
335impl<U> MulAssign<TensorRank0> for Quantity<U> {
336    fn mul_assign(&mut self, tensor_rank_0: TensorRank0) {
337        self.0 *= tensor_rank_0
338    }
339}
340
341impl<U> Div<TensorRank0> for Quantity<U> {
342    type Output = Self;
343    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
344        Self::new(self.0 / tensor_rank_0)
345    }
346}
347
348impl<U> DivAssign<TensorRank0> for Quantity<U> {
349    fn div_assign(&mut self, tensor_rank_0: TensorRank0) {
350        self.0 /= tensor_rank_0
351    }
352}
353
354impl<U> Mul<Quantity<U>> for TensorRank0 {
355    type Output = Quantity<U>;
356    fn mul(self, quantity: Quantity<U>) -> Self::Output {
357        Quantity::new(self * quantity.0)
358    }
359}
360
361/// Scaling a bare scalar by a dimensionless quantity leaves it bare.
362impl Mul<Quantity<Dimensionless>> for &TensorRank0 {
363    type Output = TensorRank0;
364    fn mul(self, quantity: Quantity<Dimensionless>) -> Self::Output {
365        self * quantity.0
366    }
367}
368
369impl<U, V> Mul<Quantity<V>> for Quantity<U>
370where
371    U: UnitMul<V>,
372{
373    type Output = Quantity<<U as UnitMul<V>>::Output>;
374    fn mul(self, quantity: Quantity<V>) -> Self::Output {
375        Quantity::new(self.0 * quantity.0)
376    }
377}
378
379impl<U, V> Mul<Quantity<V>> for &Quantity<U>
380where
381    U: UnitMul<V>,
382{
383    type Output = Quantity<<U as UnitMul<V>>::Output>;
384    fn mul(self, quantity: Quantity<V>) -> Self::Output {
385        Quantity::new(self.0 * quantity.0)
386    }
387}
388
389impl<U, V> Div<Quantity<V>> for &Quantity<U>
390where
391    U: UnitDiv<V>,
392{
393    type Output = Quantity<<U as UnitDiv<V>>::Output>;
394    fn div(self, quantity: Quantity<V>) -> Self::Output {
395        Quantity::new(self.0 / quantity.value())
396    }
397}
398
399impl<U, V> Div<Quantity<V>> for Quantity<U>
400where
401    U: UnitDiv<V>,
402{
403    type Output = Quantity<<U as UnitDiv<V>>::Output>;
404    fn div(self, quantity: Quantity<V>) -> Self::Output {
405        Quantity::new(self.0 / quantity.0)
406    }
407}
408
409impl<U, V> Mul<&Quantity<V>> for Quantity<U>
410where
411    U: UnitMul<V>,
412{
413    type Output = Quantity<<U as UnitMul<V>>::Output>;
414    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
415        self * *quantity
416    }
417}
418
419impl<U, V> Mul<&Quantity<V>> for &Quantity<U>
420where
421    U: UnitMul<V>,
422{
423    type Output = Quantity<<U as UnitMul<V>>::Output>;
424    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
425        *self * *quantity
426    }
427}
428
429impl<V> Mul<&Quantity<V>> for TensorRank0 {
430    type Output = Quantity<V>;
431    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
432        Quantity::new(self * quantity.0)
433    }
434}
435
436impl<V> Mul<&Quantity<V>> for &TensorRank0 {
437    type Output = Quantity<V>;
438    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
439        Quantity::new(self * quantity.0)
440    }
441}
442
443impl Add<TensorRank0> for Quantity<Dimensionless> {
444    type Output = Self;
445    fn add(self, tensor_rank_0: TensorRank0) -> Self::Output {
446        Self::new(self.0 + tensor_rank_0)
447    }
448}
449
450impl Add<Quantity<Dimensionless>> for TensorRank0 {
451    type Output = Quantity<Dimensionless>;
452    fn add(self, quantity: Quantity<Dimensionless>) -> Self::Output {
453        Quantity::new(self + quantity.0)
454    }
455}
456
457impl Sub<TensorRank0> for Quantity<Dimensionless> {
458    type Output = Self;
459    fn sub(self, tensor_rank_0: TensorRank0) -> Self::Output {
460        Self::new(self.0 - tensor_rank_0)
461    }
462}
463
464impl Sub<Quantity<Dimensionless>> for TensorRank0 {
465    type Output = Quantity<Dimensionless>;
466    fn sub(self, quantity: Quantity<Dimensionless>) -> Self::Output {
467        Quantity::new(self - quantity.0)
468    }
469}
470
471impl PartialEq<TensorRank0> for Quantity<Dimensionless> {
472    fn eq(&self, tensor_rank_0: &TensorRank0) -> bool {
473        &self.0 == tensor_rank_0
474    }
475}
476
477impl PartialEq<Quantity<Dimensionless>> for TensorRank0 {
478    fn eq(&self, quantity: &Quantity<Dimensionless>) -> bool {
479        self == &quantity.0
480    }
481}
482
483impl PartialOrd<TensorRank0> for Quantity<Dimensionless> {
484    fn partial_cmp(&self, tensor_rank_0: &TensorRank0) -> Option<std::cmp::Ordering> {
485        self.0.partial_cmp(tensor_rank_0)
486    }
487}
488
489impl PartialOrd<Quantity<Dimensionless>> for TensorRank0 {
490    fn partial_cmp(&self, quantity: &Quantity<Dimensionless>) -> Option<std::cmp::Ordering> {
491        self.partial_cmp(&quantity.0)
492    }
493}
494
495impl<U> Div<Quantity<U>> for TensorRank0
496where
497    U: UnitInv,
498{
499    type Output = Quantity<<U as UnitInv>::Output>;
500    fn div(self, quantity: Quantity<U>) -> Self::Output {
501        Quantity::new(self / quantity.0)
502    }
503}
504
505impl<U> std::iter::Sum for Quantity<U> {
506    fn sum<I>(iter: I) -> Self
507    where
508        I: Iterator<Item = Self>,
509    {
510        Self::new(iter.map(|quantity| quantity.0).sum())
511    }
512}
513
514impl<'a, U> std::iter::Sum<&'a Quantity<U>> for Quantity<U> {
515    fn sum<I>(iter: I) -> Self
516    where
517        I: Iterator<Item = &'a Quantity<U>>,
518    {
519        Self::new(iter.map(|quantity| quantity.0).sum())
520    }
521}
522
523impl<U> Erase for Quantity<U> {
524    type Erased = TensorRank0;
525    fn erase(&self) -> &Self::Erased {
526        &self.0
527    }
528}
529
530impl<U> Tensor for Quantity<U> {
531    type Item = Self;
532    type Unit = U;
533    fn error_count_zero(&self, tol_abs: TensorRank0, tol_rel: TensorRank0) -> Option<usize> {
534        self.0.error_count_zero(tol_abs, tol_rel)
535    }
536    fn error_count(
537        &self,
538        other: &Self,
539        tol_abs: TensorRank0,
540        tol_rel: TensorRank0,
541    ) -> Option<usize> {
542        self.0.error_count(&other.0, tol_abs, tol_rel)
543    }
544    fn full_contraction(&self, quantity: &Self) -> TensorRank0 {
545        self.0 * quantity.0
546    }
547    fn is_zero(&self) -> bool {
548        self.0 == 0.0
549    }
550    fn iter(&self) -> impl Iterator<Item = &Self::Item> {
551        std::slice::from_ref(self).iter()
552    }
553    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item> {
554        std::slice::from_mut(self).iter_mut()
555    }
556    fn len(&self) -> usize {
557        1
558    }
559    fn norm_inf(&self) -> Quantity<U> {
560        Self::new(self.0.abs())
561    }
562    fn norm_l1(&self) -> Quantity<U> {
563        Self::new(self.0.abs())
564    }
565    fn norm_p_sum(&self, p: TensorRank0) -> TensorRank0 {
566        self.0.abs().powf(p)
567    }
568    fn size(&self) -> usize {
569        1
570    }
571    fn sub_abs(&self, other: &Self) -> Self {
572        Self::new((self.0 - other.0).abs())
573    }
574    fn sub_rel(&self, other: &Self) -> Self {
575        Self::new(self.0.sub_rel(&other.0))
576    }
577}
578
579impl<U> TensorArray for Quantity<U> {
580    type Array = TensorRank0;
581    type Item = Self;
582    fn as_array(&self) -> Self::Array {
583        self.0
584    }
585    fn identity() -> Self {
586        Self::new(1.0)
587    }
588    fn zero() -> Self {
589        Self::new(0.0)
590    }
591}
592
593impl<U> FiniteDifference for Quantity<U> {
594    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
595        self.0.error_fd(&comparator.0, epsilon)
596    }
597}
598
599impl<U, const N: usize> FiniteDifference for TensorList<Quantity<U>, N> {
600    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
601        error_fd_over(self.iter().zip(comparator.iter()), epsilon)
602    }
603}
604
605impl<U, const M: usize, const N: usize> FiniteDifference
606    for TensorList<TensorList<Quantity<U>, N>, M>
607{
608    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
609        error_fd_over(
610            self.iter()
611                .zip(comparator.iter())
612                .flat_map(|(entry, comparator_entry)| entry.iter().zip(comparator_entry.iter())),
613            epsilon,
614        )
615    }
616}
617
618fn error_fd_over<'a, U: 'a>(
619    entries: impl Iterator<Item = (&'a Quantity<U>, &'a Quantity<U>)>,
620    epsilon: TensorRank0,
621) -> Option<(bool, usize)> {
622    let error_count = entries
623        .filter_map(|(entry, comparator_entry)| entry.error_fd(comparator_entry, epsilon))
624        .map(|(_, count)| count)
625        .sum();
626    if error_count > 0 {
627        Some((true, error_count))
628    } else {
629        None
630    }
631}
632
633impl<U> Solution for Quantity<U> {
634    fn decrement_from(&mut self, _other: &Vector) {
635        unimplemented!()
636    }
637    fn decrement_from_chained(&mut self, _other: &mut Vector, _vector: &Vector) {
638        unimplemented!()
639    }
640}
641
642impl<U> Hessian for Quantity<U> {
643    fn quadratic_form(&self, vector: &Vector) -> TensorRank0 {
644        self.0 * vector[0] * vector[0]
645    }
646    fn entry(&self, _row: usize, _column: usize) -> TensorRank0 {
647        unimplemented!()
648    }
649    fn fill_into(self, _square_matrix: &mut SquareMatrix) {
650        unimplemented!()
651    }
652}
653
654impl<U> Jacobian for Quantity<U> {
655    fn fill_into(&self, _vector: &mut Vector) {
656        unimplemented!()
657    }
658    fn fill_into_chained(self, _other: Vector, _vector: &mut Vector) {
659        unimplemented!()
660    }
661}
662
663impl<U> Sub<Vector> for Quantity<U> {
664    type Output = Self;
665    fn sub(self, _vector: Vector) -> Self::Output {
666        unimplemented!()
667    }
668}
669
670impl<U> Sub<&Vector> for Quantity<U> {
671    type Output = Self;
672    fn sub(self, _vector: &Vector) -> Self::Output {
673        unimplemented!()
674    }
675}
676
677impl<U> From<Vector> for Quantity<U> {
678    fn from(_vector: Vector) -> Self {
679        unimplemented!()
680    }
681}
682
683impl<U, T> Differentiate<T> for Quantity<U>
684where
685    U: UnitDiv<T>,
686{
687    type Derivative = Quantity<<U as UnitDiv<T>>::Output>;
688}
conspire/math/tensor/quantity/mod.rs

conspire/math/tensor/quantity/
mod.rs
